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ABSTRACT 


This  research  addresses  the  problem  of  acquiring  a  time  series  of  magnetic  reso¬ 
nance  images  with  both  high  spatial  and  temporal  resolutions.  Specifically,  we  sys¬ 
tematically  investigate  the  advantages  and  limitations  of  reduced-encoding  imaging 
using  a  priori  constraints.  This  study  reveals  that  if  the  available  a  priori  information 
is  a  reference  image,  direct  use  of  this  information  to  “optimize”  data  acquisition  using 
the  existing  wavelet  transform  or  singular  value  decomposition  schemes  can  under¬ 
mine  the  capability  to  detect  new  image  features.  However,  proper  incorporation  of 
the  a  priori  information  in  the  image  reconstruction  step  can  significantly  reduce  the 
resolution  loss  associated  with  reduced-encoding.  For  Fourier-encoded  data,  we  have 
shown  that  the  generalized-series  (GS)  model  is  an  effective  mathematical  framework 
for  carrying  out  the  constrained  reconstruction  step. 

Several  techniques  are  proposed  in  this  dissertation  to  improve  the  basis  functions 
of  the  GS  model  by  introducing  dynamic  information.  The  two-reference  reduced- 
encoding  imaging  by  generalized-series  reconstruction  (TRIGR)  method  suppresses 
background  information  through  the  use  of  a  second  high-resolution  reference  image. 
A  second  technique  injects  information  from  the  dynamic  data  into  the  GS  basis 
functions,  as  opposed  to  deriving  them  solely  from  the  reference  information.  These 
techniques  allow  the  GS  basis  functions  to  more  accurately  represent  the  areas  of 
dynamic  change.  Finally,  motion  that  occurs  between  the  acquisition  of  the  reference 
and  dynamic  data  sets  can  render  the  reference  information  useless  as  a  constraint 
for  image  reconstruction.  A  motion  compensation  method  is  proposed  which  uses  a 
similarity  norm  to  accurately  detect  the  motion  in  spite  of  contrast  changes  and  the 
low-resolution  nature  of  the  dynamic  data. 
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CHAPTER  1 


INTRODUCTION 

1.1  Problem  Statement 

This  research  is  about  dynamic  magnetic  resonance  imaging  (MRI).  Although 
dynamic  imaging  may  mean  different  things  in  different  applications,  from  the  data 
acquisition  standpoint,  it  can  always  be  characterized  as  the  collection  of  a  sequence 
of  images  /i(r),  J2(r)  . . .  /„(r).  Sometimes,  this  type  of  experiment  is  termed  time- 
sequential  imaging.  The  interimage  variations  AI  =  /„  —  are  commonly  called 
dynamic  changes.  For  this  study,  we  assume  that  the  dynamic  changes  may  or  may  not 
be  time-dependent  per  se.  For  example,  in  diffusion-weighted  imaging,  the  dynamic 
changes  do  not  reflect  an  underlying  time  variation  in  the  object,  but  are  due  to 
manipulation  of  the  data  acquisition  procedure  as  the  image  sequence  progresses.  On 
the  other  hand,  the  dynamic  changes  can  relate  to  a  time-dependent  change  in  the 
object  being  imaged,  such  as  that  which  occurs  when  using  MRI  to  guide  the  insertion 
of  a  biopsy  needle. 

Dynamic  MRI  is  becoming  an  increasingly  important  area  of  research  due  to  the 
many  practical  applications.  A  particular  example  of  interest  is  MR  mammography  in 
which  a  time-series  of  images  of  the  breast  is  taken  to  monitor  the  wash-in/wash-out 
of  an  injected  contrast  agent.  The  interest  in  this  lies  in  the  possibility  that  the  tem¬ 
poral  shape  of  the  enhancement  curve,  as  well  as  the  spatial  pattern  of  enhancement 
in  a  lesion,  can  determine  noninvasively  whether  a  lesion  is  benign  or  malignant  [1-9]. 
To  capitalize  on  the  period  of  greatest  differentiation  between  malignant  and  benign 
lesions,  a  sequence  of  images  must  be  acquired  during  the  first  1  or  2  minutes  follow¬ 
ing  contrast  injection  [10,11],  leading  to  the  requirement  of  high  temporal  resolution. 
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In  addition,  high  spatial  resolution  in  three  dimensions  is  imperative  to  allow  the  vi¬ 
sualization  of  very  small  tumors  with  full  coverage  of  the  breast.  A  dynamic  imaging 
method  that  could  deliver  simultaneously  high  temporal  and  spatial  resolutions  could 
detect  cancerous  lesions  at  an  earlier  stage,  thus  improving  the  patient’s  prognosis. 
Additionally,  a  method  of  noninvasive  lesion  characterization  could  reduce  the  physi¬ 
cal  and  mental  toll  on  the  patient,  as  well  as  the  financial  cost.  The  ability  to  acquire 
rapid,  high-quality  images  would  also  have  application  in  contrast-enhanced  imaging 
of  other  types  of  cancers  for  detection  [12, 13],  monitoring  the  eflPects  of  treatment  [14] 
and  watching  for  recurrence  [12]. 

The  focus  of  this  research  is  obtaining  high  spatial  and  temporal  resolutions  simul¬ 
taneously.  Specifically,  we  investigate  how  to  obtain  high-resolution  image  sequences 
with  a  reduced  number  of  dynamic  data.  Practical  issues  related  to  motion,  signal- 
to-noise  ratio  and  resolution  capability  are  also  addressed  in  this  dissertation. 

1.2  Background 
1.2.1  Signal  expression 

The  signal  activated  from  a  sample  (or  a  spin  system)  after  a  pulse  excitation  can 
be  written  in  the  following  form  [15]: 

d{t)  =  r  (1.1) 

J  —  00 

where  d{t)  is  the  activated  signal,  W{f)  is  a  window  function  that  represents  the 
excitation  sensitivity,  M{f)  is  the  spatial  distribution  of  spins  in  the  object,  7  is  a 
proportionality  constant  and  a;(f}  is  the  frequency  of  the  activated  signal.  Typically, 
W (f)  defines  a  planar  slice  of  the  object  [16],  but  more  sophisticated  spatially  selective 
excitation  schemes  are  possible  such  as  those  used  in  non-Fourier  encoding  methods. 
Note  that  the  received  signal  in  Eq.  (1.1)  is  the  integral  over  the  entire  excited  region 
of  the  object. 
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1.2.2  The  issue  of  imaging  time 

In  order  to  create  an  image  from  the  activated  signal,  it  is  necessary  to  encode 
spatial  information  into  the  signal  [17].  Frequency-encoding  makes  the  frequency  of 
the  acquired  signal  depend  on  the  location  of  the  spins  as 

Aa;(r/e)  =  7<5/er/e,  (1-2) 

where  Au{rfe)  is  the  change  in  frequency  as  a  function  of  position  along  the  frequency¬ 
encoding  direction  and  G/e  is  the  applied  frequency-encoding  gradient.  The  phase¬ 
encoding  method  encodes  the  spatial  information  into  the  initial  phase  of  the  received 
signal  as 

~  'yG p^Tp^Tp^,  (1-3) 

where  (p^Vpe)  is  the  phase  as  a  function  of  position  along  the  phase-encoding  direction, 
Gpe  is  the  applied  phase-encoding  gradient  and  Tp^  is  the  phase-encoding  time  period. 
The  encoding  process  can  be  described  using  the  popular  A:-space  notation  [18]  in 
which  the  A:-space  variable  k  is  defined  as 

kt) = (1.4) 

where  G  is  the  applied  gradient  magnetic  field.  Using  this  notation,  the  imaging 
equation  is  given  by  the  Fourier  transform  as 

d{k)  =  /  (1.5) 

7—00 

where  d{k)  is  the  acquired  signal  and  /(f)  is  the  desired  image  function. 

To  create  a  good  quality  image,  measurements  must  be  made  from  many  locations 
in  A:-space.  In  the  conventional  Fourier  imaging  technique,  a  combination  of  the 
phase  and  frequency-encoding  methods  discussed  above  is  used.  During  each  data 
acquisition  period,  a  set  of  data  points  is  acquired  using  frequency  encoding  as 

f^fein)  =  ^GfeuAt  — ^  <  n  <  (1.6) 
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Figure  1.1:  Coverage  of  fc-Space:  (a)  standard  Fourier  encoding,  (b)  pro¬ 
jection  reconstruction,  (c)  echo  planar  and  fast  spin-echo,  (d)  spiral  scan¬ 
ning  and  (e)  concentric  scanning. 

where  Nje  is  the  number  of  points  acquired  in  the  frequency-encoding  direction  and 
At  is  the  time  increment  between  samples.  For  the  phase-encoding  direction,  the  value 
of  the  phase- encoding  gradient  is  incremented  between  each  of  the  N^,  encodings  as 

kpe{n)  =  ^nAGpeFpe  <n<^,  (1.7) 

where  AGpe  is  the  step  in  phase-encoding  gradient.  The  resulting  coverage  of  /j-space 
is  shown  in  Fig.  1.1(a).  Another  example  is  projection  reconstruction  imaging,  which 
uses  frequency  encoding  for  two  dimensions.  In  this  case,  A:-space  is  covered  in  a  radial 
fashion  as  shown  in  Fig.  1.1(b)  by  incrementing  two  frequency-encoding  gradients  as 

kfei{n)  =  cos{nA(f))t  and  (1.8a) 

kfe2{n)  =  ^Gsin(nA(/))t  <  n  <  ^,  (1.8b) 

where  G  is  a  constant  and  A(j>  is  the  angular  increment  in  A:-space. 

No  matter  what  coverage  of  A:-space  is  selected,  the  time  to  acquire  an  image 
depends  upon  the  number  of  encodings  that  are  required.  If  Tr  is  the  repetition  time 
or  the  time  for  one  encoding  and  Ne  excitations  are  necessary,  the  total  imaging  time 
will  be 

T  =  N^Tr.  (1.9) 

In  many  cases,  the  minimum  value  of  Tr  is  limited  by  the  Ti  tissue  relaxation  param¬ 
eter  in  which  case,  an  increase  in  temporal  resolution  would  have  to  come  through  a 
reduction  in  the  number  of  acquired  encodings. 
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1.2.3  The  issue  of  image  resolution 

Describing  the  MR  imaging  process  in  terms  of  linear  system  theory,  the  resulting 
image  /(f)  will  be  related  to  the  original  object  /(f)  through  a  convolution  as 


/(f)  =  /(f)  *  h{f), 


(1.10) 


where  h{f)  is  the  point  spread  function  (PSF)  of  the  imaging  process.  The  width  of 
the  PSF  defines  the  separation  required  between  two  points  in  the  object  such  that 
they  can  be  resolved  in  the  image  and  thus,  gives  a  measure  of  the  spatial  resolution 
of  the  image.  The  width  of  the  PSF  is  often  described  in  terms  of  the  width  w  of  an 
equivalent  rectangle  as 

J  /i(0)rect  ^  dr  =  h{r)dr,  (l-H) 


where 


rect(r)  = 


r  < 


0  else. 

For  example,  if  the  standard  Fourier  imaging  method  is  used,  the  PSF  will  be 

Mr)  = 


(1.12) 


(1.13) 


sin(7rAA:r) 

where  Ng  is  the  number  of  A;-space  data  points  and  Ak  is  the  step  size  in  fc-space. 
Therefore,  the  spatial  resolution  of  the  resulting  image  is 

1 


Ar  = 


NgAk’ 


(1.14) 


where  Ar  is  the  pixel  size  in  the  image.  By  comparing  Eqs.  (1.9)  and  (1.14),  it  can 
be  seen  that  an  improvement  in  the  temporal  resolution  through  a  reduced  number 
of  phase  encodings  will  be  accompanied  by  a  commensurate  loss  in  spatial  resolution 
for  the  conventional  Fourier  imaging  technique. 


1.2.4  Fast  imaging 

For  many  dynamic  imaging  applications,  fast  imaging  techniques  are  necessary  to 
provide  adequate  temporal  resolution.  It  is  clear  from  Eq.  (1.9)  that  one  can  increase 
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the  temporal  resolution  of  a  sequence  of  images  by  either  shortening  Tr  or  reducing 
iVg.  Techniques  that  use  the  former  strategy  are  called  fast-scan  techniques  and  the 
latter  are  called  reduced-scan  techniques. 

Fast-scan  imaging  techniques  try  to  acquire  a  full  set  of  data  in  a  time  that  is  short 
compared  to  the  dynamic  process.  The  /j-space  coverage  of  several  of  these  techniques 
is  illustrated  in  Fig.  l.l(c)-(e).  Echo  planar  [19,20]  uses  a  rapidly  switching  frequency 
encoding  gradient  to  cover  A:-space  in  a  rectilinear  fashion  during  data  acquisition.  The 
fast  spin-echo  technique  [21-23]  obtains  the  same  /c-space  coverage  as  echo  planar,  but 
uses  multiple  refocussing  pulses  instead  of  the  gradient  switching.  Because  the  data 
for  these  methods  are  acquired  on  a  rectilinear  grid,  image  reconstruction  can  be  easily 
performed  using  the  fast  Fourier  transform  (FFT).  Spiral  scanning  [24-26]  applies 
time- varying  frequency  encoding  gradients  during  data  acquisition  to  traverse  a  spiral 
trajectory  in  A:-space.  Image  reconstruction  for  spiral  scanning  requires  interpolation 
to  a  Cartesian  grid  followed  by  the  FFT  [27,28].  With  this  reconstruction  method, 
it  is  necessary  to  know  the  exact  A:-space  location  of  the  data  points  to  avoid  image 
artifacts,  often  necessitating  the  measurement  of  the  actual  fc-space  trajectory  due 
to  imperfect  gradient  performance  [29,30].  Concentric  scanning  [31,32]  also  uses 
a  modulated  frequency  encoding  gradient,  but  acquires  circles  of  data  in  A;-space 
instead  of  a  spiral  of  data.  For  all  of  these  techniques,  an  image  can  be  generated 
with  a  single  excitation  or  interleaved  segments  of  the  required  data  can  be  acquired 
wdth  a  few  excitations.  Although  fast-scan  techniques  are  capable  of  acquiring  an 
image  in  a  short  time,  they  may  require  specialized  hardware  for  implementation. 
In  addition,  they  may  have  reduced  contrast  between  tissues  and  a  reduced  signal- 
to-noise  ratio  (SNR).  For  dynamic  imaging  applications,  there  may  also  be  power 
deposition  problems  due  to  the  multiple  images  that  are  acquired. 

In  contrast  to  the  fast-scan  methods,  the  reduced-scan  methods  try  to  improve 
the  temporal  resolution  by  reducing  Ng.  These  techniques  are  motivated  by  the 
observation  that  in  many  dynamic  imaging  applications,  the  dynamic  process  induces 
a  contrast  modulation  on  the  underlying  static  high-resolution  morphology.  They 
exploit  this  fact  to  reduce  the  redundancy  in  the  data  acquisition  process  to  improve 
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Reference  Dynamic  Imaging  Period 

Figure  1.2:  Reduced- Encoding  Data  Acquisition:  A  high-resolution  refer¬ 
ence  data  set  is  followed  by  a  series  of  reduced-encoding  data  sets  during 
the  dynamic  imaging  period. 

the  temporal  resolution.  The  reduced-scan  data  acquisition  strategy  is  illustrated 
in  Fig.  1.2.  A  high-resolution  reference  data  set  is  acquired,  followed  by  a  sequence 
of  reduced-encoding  dynamic  data  sets.  If  M  and  N  encodings  are  acquired  for  the 
reference  and  dynamic  data  sets,  respectively,  the  improvement  in  temporal  resolution 
for  the  dynamic  images  will  be  M/N.  However,  as  discussed  before,  this  will  result  in 
a  loss  of  spatial  resolution  with  the  conventional  Fourier  reconstruction  method.  To 
avoid  this  loss  of  spatial  resolution,  many  of  the  reduced-encoding  dynamic  imaging 
methods  use  a  priori  information  from  the  reference  image  at  some  point  in  the 
imaging  process  to  reduce  the  number  of  encodings  required  per  dynamic  image. 

Several  of  the  reduced-scan  techniques,  such  as  wavelet  encoding  methods  [33-57] 
and  singular  value  decomposition  methods  [58-70],  use  the  a  priori  information  during 
the  data  acquisition  process  to  determine  a  more  optimal  truncated  basis  set  than 
the  infinite  complex  exponentials  of  the  Fourier  encoding  set.  The  hope  is  that  the 
non-Fourier  basis  vectors  will  be  better  able  to  encode  the  information  in  the  object 
with  a  small  number  of  basis  functions.  Other  reduced-scan  techniques  use  the  a 
priori  information  as  a  constraint  for  extrapolation  of  the  missing  dynamic  data 
during  image  reconstruction  [71].  These  include  the  reduced-encoding  imaging  by 
generalized-series  reconstruction  (RIGR)  method  [72,73]  and  the  keyhole  or  data- 
replacement  technique  (DRT)  [74,75].  The  reduced-scan  methods  will  be  discussed 
in  more  detail  in  Chapter  2. 
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1.3  Summary  of  Results 


The  key  issue  addressed  in  this  research  is  how  to  obtain  simultaneously  high  spa¬ 
tial  and  temporal  resolutions  within  the  reduced-encoding  dynamic  imaging  frame¬ 
work.  To  avoid  the  loss  of  spatial  resolution  that  occurs  with  reduced-encoding  Fourier 
imaging,  many  methods  make  use  of  a  priori  information  to  reduce  the  number  of 
encodings  that  are  necessary  for  the  dynamic  images.  In  this  dissertation,  it  was 
determined  that  if  the  a  priori  information  that  is  available  is  a  reference  image,  di¬ 
rect  use  of  this  information  to  “optimize”  data  acquisition  using  the  existing  wavelet 
transform  or  singular  value  decomposition  schemes  can  undermine  the  capability  to 
detect  new  image  features.  Incorporation  of  the  a  priori  information  in  the  image 
reconstruction  step  is  preferable.  We  have  also  shown  that  the  generalized  series  (GS) 
model  of  the  RIGR  technique  is  a  better  way  to  combine  the  reference  and  dynamic 
data  for  image  reconstruction  than  the  keyhole  method. 

Several  techniques  are  proposed  in  this  dissertation  to  improve  the  GS  basis  func¬ 
tions  by  introducing  dynamic  information.  The  two-reference  reduced-encoding  imag¬ 
ing  by  generalized-series  reconstruction  (TRIGR)  method  suppresses  background  in¬ 
formation  through  the  use  of  a  second  high-resolution  reference  image.  This  allows 
the  GS  basis  functions  to  more  accurately  represent  the  areas  of  dynamic  change. 
A  second  technique  injects  information  from  the  dynamic  data  into  the  GS  basis 
functions  as  opposed  to  deriving  them  solely  from  the  reference  information.  The 
resulting  GS  basis  functions  more  closely  resemble  those  that  would  be  derived  from 
the  dynamic  image  itself.  Finally,  motion  that  occurs  between  the  acquisition  of  the 
reference  data  set  and  the  dynamic  data  sets  can  render  the  reference  information 
useless  as  a  constraint  for  image  reconstruction.  This  is  a  difficult  problem  to  address, 
because  the  dynamic  changes  may  alter  the  appearance  of  the  image  significantly,  pos¬ 
ing  problems  for  both  navigator-based  techniques  and  registration  algorithms.  The 
proposed  method  uses  a  similarity  norm  to  accurately  detect  the  motion,  in  spite  of 
the  contrast  changes,  and  can  significantly  reduce  motion  artifacts  in  GS  images. 
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1.4  Organization  of  the  Dissertation 


The  remainder  of  the  dissertation  is  organized  as  follows.  Chapters  2  and  3  investi¬ 
gate  issues  involved  in  data  acquisition  and  image  reconstruction  for  reduced-encoding 
dynamic  imaging.  Chapter  4  formulates  a  method  to  use  the  additional  information 
in  a  second  high-resolution  reference  image  to  suppress  the  background  information  in 
the  GS  basis  functions.  Chapter  5  details  a  method  to  use  the  edge  information  from 
the  reference  image  with  the  contrast  information  from  the  dynamic  data  to  improve 
the  resulting  dynamic  images.  Chapter  6  introduces  a  motion  detection  scheme  for 
reduced-encoding  dynamic  imaging  applications.  Chapter  7  contains  suggestions  for 
future  work  and  the  conclusions  drawn  from  this  research. 
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CHAPTER  2 


DATA  ACQUISITION 


For  convenience,  we  shall  denote  the  MRI  data  in  terms  of  an  inner  product  as 

d(n)  =  (/,  e(n))  n  €  A/",  (2.1) 

where  d{n)  is  the  acquired  data,  I  is  the  object  function  and  e(n)  are  the  encoding 
vectors.  For  reduced-encoding  dynamic  imaging,  the  key  issue  is  the  selection  of  a 
truncated  set  of  basis  functions  e(n),  n  G  Mdyn  C  Af,  to  best  represent  the  dynamic 
image  I.  With  the  typical  Fourier  reduced-encoding  scheme,  a  priori  information  is 
not  used  to  select  the  truncated  set  of  encoding  vectors.  If  it  is  desired  to  inject  a 
priori  information  into  the  data  acquisition,  two  current  methods  for  doing  that  are 
the  wavelet  and  singular  value  decomposition  (SVD)  encoding  methods.  The  wavelet 
encoding  scheme  uses  the  a  priori  information  only  to  guide  the  truncation;  whereas, 
the  SVD  method  uses  the  information  to  both  select  and  truncate  the  encoding  vector 
set.  For  this  research,  the  wavelet  and  SVD  encoding  methods  were  investigated  in 
comparison  to  Fourier  encoding  for  reduced-encoding  dynamic  imaging  applications. 

2.1  Fourier  Encoding 

With  the  Fourier  encoding  method,  the  encoding  vectors  are  the  Fourier  basis 
set  of  infinite  complex  exponentials  where  k  is  a  spatial  frequency  variable. 

Although  it  would  be  possible  to  adaptively  truncate  the  Fourier  encoding  vectors  that 
are  acquired  for  an  object,  typically  no  a  priori  information  is  used  in  the  Fourier- 
encoding  data  acquisition  process.  Given  that  no  a  priori  information  is  available 
about  the  object,  it  is  optimal  to  sample  data  from  the  center  of  ^-space  because  of 
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the  well-known  fact  that  /c-space  data  decays  as  \/k  for  practical  image  functions. 
This  scheme  can  also  be  motivated  by  the  observation  that  the  central  A;-space  data 
contributes  the  bulk  of  the  contrast  information  to  the  resulting  image;  whereas,  data 
in  the  outer  region  of  /:-space  largely  represents  the  edge  information.  Therefore,  with 
Fourier  encoding,  the  reduced-encoding  data  set  that  is  collected  is  described  by 

d{nLk)  =  l°°  ^  <  n  <  y  -  1.  (2.2) 

2.2  Wavelet  Encoding 

In  contrast  to  the  infinite  complex  exponentials  of  the  Fourier  basis  set,  the  wavelet 
basis  functions  are  localized  in  both  the  frequency  and  spatial  domains.  The  basis  set 
is  formed  from  dilations  and  translations  of  a  scaling  function  4>,  which  extracts  an 
approximation  of  the  object  function,  and  the  associated  wavelet  </?,  which  extracts 
the  details  of  the  object  function.  To  form  a  wavelet  basis  set,  the  ^  must  be  con¬ 
tinuously  differentiable  and  must  satisfy  I^(r)|  =  0{r~'^)  and  |</''(r)|  =  0(r~^)  [76]. 
Alternatively,  the  wavelet  function  must  satisfy  the  admissibility  condition  [77] 

J  1^1”^  (2-3) 

The  commonly  used  wavelets  can  be  classified  into  two  general  categories:  orthogonal 
and  biorthogonal.  An  orthogonal  wavelet  transform  requires  two  basis  sets.  For  a 
J-level  discrete  dyadic  orthogonal  wavelet  transform,  the  set  of  these  functions  can 
be  expressed  as 

(f>-j,k{r)  =  4>{2-^r-k)  A:  =  1, 2 . . . ,  M2-^  and  (2.4a) 

<^,y(r)  =  V^ip{2^r-l)  j-  =  -l,-2...,-J;/  =  l,2...,M2^  (2.4b) 

where  ^  and  (/?  are  the  scaling  function  and  wavelet,  respectively,  and  M  is  the 
number  of  points  across  the  object  function.  The  other  category  is  biorthogonal 
wavelets,  which  have  the  advantage  of  linear  phase.  (It  has  been  shown  that  tradi¬ 
tional  wavelets  can  only  be  orthogonal  and  linear  phase  simultaneously  for  the  Haar 
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wavelet  basis  [77].)  For  biorthogonal  wavelets,  four  basis  sets  are  required  because 
the  analysis  and  synthesis  scaling  functions  (wavelets)  cannot  be  the  same  function  as 
in  the  orthogonal  case.  For  a  J-level  discrete  dyadic  biorthogonal  wavelet  transform, 
the  set  of  these  functions  can  be  expressed  as- 


=  \/2  ^  ^[2  —  kij 

k  =  l,2...,M2-’, 

(2.5a) 

=  p  (2^r  — 

i  =  -i,-2.. 

.,-J;/  =  1,2., 

..,M2^ 

(2.5b) 

=  ^  (2-^r  -  k) 

A:  =  1, 2 . . . ,  M2--^ 

and 

(2.5c) 

=  <p  (2^r  -  l) 

.,-J;/  =  1,2.. 

.,M2\ 

(2.5d) 

where  (j)  and  v?  of  Eq.  (2.4)  are  the  analysis  scaling  function  and  wavelet,  respectively, 
and  ^  and  ip  are  the  synthesis  scaling  function  and  wavelet,  respectively.  Note  that 
the  basis  functions  are  down-sampled  by  a  factor  of  two  at  each  level  due  to  the 
dilation,  and  hence,  the  transform  will  not  have  redundant  information.  Therefore,  a 
signal  of  M  samples  can  be  represented  by  a  J-level  discrete  wavelet  transform  of  M 
total  coefficients  scaling  coefficients  and  MYljd-i  2^  wavelet  coefficients). 

By  wavelet  theory,  the  wavelet  coefficients  that  are  the  largest  are  the  most  im¬ 
portant,  and  these  tend  to  be  sparse,  meaning  that  most  object  functions  can  be  well 
represented  with  a  truncated  set  of  wavelet  coefficients.  To  try  to  exploit  this  prop¬ 
erty  for  dynamic  MRI  applications,  wavelet  encoding  is  applied  in  the  following  way. 
First,  a  high-resolution  reference  image  is  collected,  typically  with  Fourier  encoding. 
This  reference  image  is  decomposed  along  the  proposed  wavelet  encoding  direction 
using  the  ID  wavelet  transform.  Based  on  a  chosen  criteria,  the  “most  important”  N 
wavelet  basis  functions  are  selected  to  be  used  as  excitation  profiles  for  the  dynamic 
imaging  period.  This  reduced  set  of  encodings  that  is  acquired  can  be  expressed  as, 
for  a  J-level  wavelet  transform, 

■^2--^{-7dyn}  —  ^Jdyn(^))2  (p  (2  '^r  —  k  =  ki,  k2 . .  ■  and  (2.6a) 

^>2.{Jdyn}  =  (/dy„(r),2V(2^r-/))  j  =  -l,-2...-J-l  =  h,l2...lN^j,{2.6h) 
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where  Ntj,  scaling  functions  and  N^j',  j  =  — 1,— 2...,— J  wavelets  are  utilized  and 

N  =  Nt  +  Ejd-i 

To  perform  wavelet  encoding  of  an  object,  one  must  explicitly  form  the  inner 
product  of  the  object  and  the  basis  functions  using  spatially  selective  RF  pulses  to 
generate  an  excitation  profile  in  the  shape  of  the  selected  basis  functions.  In  other 
words,  the  acquired  signal  is  received  from  a  scaling  function  or  wavelet-shaped  region 
of  the  object.  The  excitation  profile  is  translated  and  dilated  to  acquire  the  set  of 
data  specified  in  Eq.  (2.6).  Translation  of  the  excitation  profile  is  accomplished  by 
either  shifting  the  center  frequency  of  the  RF  pulse  [40]  or  ramping  the  phase  during 
the  pulse  [54].  Dilation  of  the  excitation  profile  can  be  achieved  by  changing  the 
gradient  strength  [40]  or  changing  the  length  of  the  i?E  pulse  [53]. 

2.3  Singular  Value  Decomposition  Encoding 

Another  non-Fourier  encoding  method  is  based  on  the  singular  value  decomposi¬ 
tion  (SVD).  The  SVD  of  a.  j  x  k  matrix  A  of  rank  r  is  defined  as 

A  =  (2.7) 

where  U  and  V  are  j  x  j  and  k  x  k  unitary  matrices,  respectively,  and  'E  is  a,  j  x  k 
diagonal  matrix.  The  columns  of  U  and  V  are  called  the  left  and  right  singular 
vectors,  respectively,  and  the  r  diagonal  entries  of  E  are  called  the  singular  values. 
The  SVD  can  also  be  written  as 

A  =  ^  ffiUivf,  (2.8) 

i=l 

where  Ui  and  Vi  are  the  columns  of  U  and  V,  respectively,  and  cJi  are  the  diagonal 
entries  of  E. 

With  the  SVD,  it  is  easy  to  choose  the  optimal  truncated  encoding  vector  set 
using  the  Eckart  and  Young  theorem  [78].  This  theorem  states  that  the  minimum 
distance  in  the  L2  norm  sense  between  a  matrix  A  of  rank  r  and  any  matrix  B  of 
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rank  less  than  or  equal  to  f  <  r  is  achieved  by  the  truncated  SVD^  which  is  defined 
as 

A  =  (2.9a) 

i=l 

=  (2.9b) 

Here  Er  is  a  diagonal  matrix  containing  the  r-most  significant  singular  values,  and  Uf 
and  Vf  are  the  matrices  consisting  of  the  left  and  right  singular  vectors,  respectively, 
that  are  associated  with  the  f-most  significant  singular  values. 

This  feature  of  the  SVD  has  been  exploited  in  many  fields  and  is  the  basis  for  its 
application  to  dynamic  MR  imaging.  The  truncated  SVD  representation  of  a  reference 
image  is  used  to  define  the  encoding  vectors  for  the  dynamic  data  acquisition,  in  the 
hopes  that  these  vectors  would  also  form  a  good  representation  for  the  dynamic 
image.  Specifically,  given  a  reference  image  /ref,  the  singular  value  decomposition  of 
the  reference  image  is  first  performed  as 

7ref  =  c/sy^.  (2.10) 

For  SVD  encoding  along  the  vertical  or  horizontal  directions,  the  Ui  or  Vi,  respectively, 
corresponding  to  the  largest  N  singular  values  are  selected  as  the  encoding  vectors 
for  the  dynamic  imaging  period.  In  the  pure  vertical  encoding  case,  for  example,  the 
dynamic  data  set  generated  for  each  dynamic  image  is  the  projection  of  the  ideal 
image  I^yn  onto  the  space  spanned  by  the  selected  left  singular  vectors.  This  can  be 
expressed  mathematically  as 

•^dyn  ==  /(dynj  (2-11) 

where  Un  is  the  matrix  constructed  from  the  N  selected  left  singular  vectors.  As 
with  wavelet  encoding,  SVD  encoding  requires  that  the  inner  product  is  explicitly 
created  using  spatially  selective  RF  pulses  to  excite  a  region  of  the  object  in  the 
shape  of  a  selected  SVD  singular  vector.  Unlike  wavelet  encoding,  which  can  use  many 
translations  and  dilations  of  a  given  excitation  profile  depending  upon  the  selected 
wavelet  basis  functions,  each  SVD  encoding  will  involve  a  different  excitation  profile. 

^Note  that  the  minimizer  of  \\A  —  B\\2  is  unique  if  and  only  if  cr^  >  <7f+i. 
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Scaling  Function  Wavelet  Function 


Figure  2.1:  Daubechies  D18  Orthogonal  Wavelet:  (a)  scaling  function  and 
(b)  wavelet 

2.4  Encoding  Method  Assessment 


The  methods  were  compared  using  computer  simulations  on  both  simulated  and 
real  MRI  data  of  dynamic  imaging  applications.  To  simulate  the  reduced  phase¬ 
encoding  data  acquisition,  A:-space  data  sets  were  generated  from  a  sequence  of  high- 
resolution  images.  A  baseline  high-resolution  data  set  was  used  as  the  reference  data, 
and  the  central  N  phase  encodings  from  the  remaining  data  sets  were  used  as  the 
dynamic  phase  encodings.  To  simulate  wavelet  encoding,  the  reference  image  was  de¬ 
composed  using  the  ID  wavelet  transform.  The  N  “most  important”  reference  wavelet 
encodings  w'ere  selected  based  on  the  squared  sum  along  the  frequency-encoding  di¬ 
rection  for  each  wavelet.  The  dynamic  data  associated  with  the  selected  wavelet 
encodings  were  generated  using  Eq.  (2.6).  To  simulate  the  SVD  encoding  method, 
the  singular  value  decomposition  of  the  reference  image  was  calculated,  and  the  left 
singular  vectors  associated  with  the  N  largest  singular  values  were  selected  as  the 
dynamic  encodings.  The  dynamic  data  were  then  generated  using  Eq.  (2.11). 

Note  that  with  the  wavelet  encoding  method,  the  choice  of  wavelet  basis  set  will 
affect  the  results  due  to  implementation  issues,  as  well  as  the  ability  to  represent  the 
dynamic  image  with  a  truncated  basis  set.  In  selecting  a  particular  wavelet  basis  set 
for  experimental  MR  encoding,  a  smoother  wavelet  is  preferred,  because  it  requires 
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Figure  2.2:  Cohen  7/9  Biorthogonal  Wavelet:  (a)  analysis  seeding  func¬ 
tion,  (b)  zuialysis  wavelet,  (c)  synthesis  scaling  function  and  (d)  synthesis 
wavelet 


a  shorter  jRF  pulse  [40]  and  reduces  the  bandwidth  requirement  of  the  pulse  [52]. 
In  addition,  a  more  accurate  wavelet-shaped  profile  can  be  excited,  which  will  result 
in  better  images  [40].  However,  there  is  a  trade-off  between  the  length  of  the  required 
RF  pulse  and  the  spatial  support  of  the  wavelet.  This  research  utilized  the  orthogonal 
Daubechies  D18  wavelet  basis  set  [77]  and  the  biorthogonal  Cohen  7/9  wavelet  basis 
set  [79],  pictured  in  Figs.  2.1  and  2.2,  respectively,  which  are  considered  by  many 
to  be  among  the  best  for  image  compression  [80].  This  dissertation  focused  on  the 
theoretical  power  of  the  basis  set  for  encoding  dynamic  images  and  chose  to  ignore 
the  long  RF  pulse  that  would  be  required  to  excite  the  profiles  associated  with  these 
wavelets. 

In  addition,  the  non-Fourier  encoding  methods  described  here  use  the  reference 
data  set  to  choose  the  encoding  vectors  for  the  entire  dynamic  imaging  period.  A 
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modification  of  this  technique  is  to  use  a  reconstructed  dynamic  image  to  select  the 
encodings  for  the  subsequent  dynamic  image.  Although  this  should  help  to  increase 
the  similarity  between  the  image  used  to  truncate  the  encoding  vector  set  and  the 
image  acquired  with  the  selected  dynamic  encodings,  it  would  reduce  the  achiev¬ 
able  temporal  resolution  due  to  the  computation  and  magnet  setup  required  between 
successive  dynamic  data  sets. 

It  is  important  to  note  that  the  simulations  in  this  study  did  not  take  into  ac¬ 
count  several  factors  that  would  degrade  the  performance  of  the  non-Fourier  encoding 
methods.  First,  it  was  assumed  that  the  spatially  selective  RF  encodings  could  be 
exactly  excited.  This  is  difficult  in  practice,  and  imperfect  excitation  will  degrade  the 
resulting  image  [40,67].  In  addition,  Ti  variations  across  the  image  will  further  distort 
the  encoded  profile  when  a  short  is  used,  which  is  especially  a  problem  when  a 
Ti  contrast  agent  is  injected.  The  use  of  RF  encoding  with  the  non-Fourier  methods 
also  limits  the  application  to  single  slice  or  3D  spin  echo  imaging,  and  it  is  not  easy 
to  implement  2D  or  thin  slab  3D  gradient  echo  sequences  [61].  Another  point  is  the 
signal-to-noise  ratio  (SNR)  loss  due  to  the  use  of  spatially  selective  excitation  for 
the  non-Fourier  encoding  methods  [40,47].  These  problems  do  not  arise  with  Fourier 
encoding,  which  uses  linear  gradients  to  encode  the  image  information  rather  than 
spatially  selective  excitation. 

As  discussed  earlier,  the  truncated  wavelet  and  SVD  representations  have  the  de¬ 
sirable  property  of  representing  an  image  well  with  fewer  encodings  than  are  necessary 
with  the  Fourier  encoding  method.  To  exploit  these  desirable  properties  for  reduced- 
encoding  dynamic  imaging,  the  current  methods  use  a  reference  image  to  guide  the 
truncation  of  the  set  of  encoding  vectors.  The  danger  in  doing  that  is  that  the  selected 
encodings  may  not  be  optimal  for  the  representation  of  the  dynamic  image  and,  at 
times,  the  reference-based  truncation  can  introduce  dangerous  artifacts.  To  illustrate 
this,  simulations  of  two  important  dynamic  imaging  applications,  a  contrast-enhanced 
dynamic  study  and  an  interventional  MRI  needle  biopsy  procedure,  are  shown  below. 

Figure  2.3  shows  the  results  of  a  contrast-enhanced  simulation  in  which  the  dy¬ 
namic  changes  involve  a  variable  rate  of  enhancement  in  each  of  the  four  “lesions,” 
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as  well  as  a  slow  overall  background  enhancement.  Figures  2.3(a)-(b)  show  the  ref¬ 
erence  and  dynamic  images,  respectively,  reconstructed  using  128  phase  encodings. 
Figures  2.3(c)-(f)  show  the  dynamic  image  reconstructed  using  16  dynamic  encodings 
using  the  Fourier,  orthogonal  wavelet,  biorthogonal  wavelet  and  SVD  encoding  meth¬ 
ods,  respectively.  The  SVD  method  results  in  a  better  delineation  of  the  boundaries 
of  the  lesions  than  the  other  methods.  However,  it  fails  to  faithfully  reproduce  the 
signal  magnitudes  in  the  lesions.  This  is  quantified  in  Fig.  2.4  which  shows  the  aver¬ 
age  signal  magnitude  in  each  of  the  lesions  as  assigned  by  the  different  methods.  Note 
that  the  SVD  methods  assign  nearly  the  same  signal  magnitude  to  all  four  lesions, 
which  is  undesirable  because  the  aim  of  this  application  is  to  accurately  track  the 
signal  magnitude  change  in  the  lesions. 

In  a  needle  biopsy  application,  the  purpose  of  MR  imaging  is  to  accurately  local¬ 
ize  the  needle  to  ensure  that  the  lesion  is  biopsied  as  opposed  to  surrounding  normal 
tissue.  The  results  of  a  needle  biopsy  simulation  are  shown  in  Fig.  2.5  in  which 
(a)-(b)  are  the  reference  and  dynamic  image,  respectively,  reconstructed  using  256 
phase  encodings.  The  dynamic  image  contains  an  additional  dark  line  feature,  sup¬ 
posedly  created  by  the  insertion  of  a  biopsy  needle.  The  dynamic  image  is  shown 
in  (c)-(f),  reconstructed  using  32  dynamic  encodings  with  the  Fourier,  orthogonal 
wavelet,  biorthogonal  wavelet  and  SVD  encoding  methods,  respectively.  Note  the 
apparent  displacement  of  the  center  of  the  needle  in  the  SVD  image  and  the  smeared 
version  of  the  needle  in  the  wavelet  reconstructions.  In  the  Fourier  image,  the  familiar 
Gibbs  ringing  results  from  the  sharp  edge  features. 

The  smearing  and  displacement  artifacts  seen  in  the  wavelet  and  SVD  images 
are  not  due  to  the  truncation  level,  but  arise  from  the  reference-based  truncation 
method.  To  illustrate  this.  Figs.  2.6(a)-(c)  were  reconstructed  using  the  optimal  32 
orthogonal  wavelet,  biorthogonal  wavelet  and  SVD  encodings  as  determined  by  the 
dynamic  image  itself.  The  greatly  improved  reconstruction  of  the  needle  over  that  in 
Figs.  2.5(d)-(f)  attest  the  nonoptimality  of  the  dynamic  encodings  selected  using  the 
reference  image  of  Fig.  2.5(a).  However,  note  the  bright  needle  ghost  artifact  in  the 
orthogonal  wavelet  case,  which  occurs  due  to  the  localized  nature  of  the  wavelet  basis 
functions. 
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Figure  2.3:  Non-Fourier  Encoding  Applied  to  Contrast-enhanced  Imaging. 
(a)-(b)  The  reference  and  dynamic  images,  respectively,  reconstructed  us¬ 
ing  128  phase  encodings.  The  remaining  images  were  reconstructed  with 
16  dynamic  encodings  using  different  encoding  methods:  (c)  Fourier,  (d) 
orthogonal  wavelet,  (e)  biorthogonal  wavelet,  and  (f)  SVD.  The  Fourier, 
wavelet  and  SVD  encoding  directions  are  vertical.  Note  that  the  SVD 
method  assigns  nearly  the  same  signal  magnitude  to  all  four  lesions. 
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Figure  2.4:  Average  Signal  Magnitude  of  the  Lesions  of  Fig.  2.3(b)-(f). 
Regions  1-4  correspond  to  the  upper-left,  upper-right,  lower-left  and  lower- 
right  lesions,  respectively,  in  Fig.  2.3. 
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Figure  2.5:  Non-Fourier  Encoding  Applied  to  Interventional  MRI:  (a)-(b) 
The  reference  and  dynamic  images,  respectively,  reconstructed  using  256 
phase  encodings.  The  remaining  images  were  reconstructed  using  32  dy¬ 
namic  encodings  with  different  encoding  techniques:  (c)  Fourier,  (d)  or¬ 
thogonal  wavelet,  (e)  biorthogonal  wavelet  and  (f)  SVD.  The  Fourier, 
wavelet  and  SVD  encoding  directions  are  horizontal.  The  arrow  indicates 
the  center  of  the  needle  track.  Note  the  apparent  displacement  of  the 
needle  center  in  the  SVD  reconstruction. 
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Figure  2.6:  Optimal  Non-Fourier  Encoding:  Dynamic  image  reconstructed 
with  32  optimal  encodings  as  derived  from  the  dynamic  image  itself  using 
(a)  orthogonal  wavelet,  (b)  biorthogonal  wavelet  and  (c)  SVD.  The  wavelet 
and  SVD  encoding  directions  are  horizontal,  and  the  arrow  indicates  the 
center  of  the  needle  track.  Note  the  improvement  over  the  images  in 
Fig.  2.5(d)-(f),  respectively,  which  were  reconstructed  using  the  32  sub- 
optimal  encodings  as  determined  from  the  reference  image  in  Fig.  2.5(a). 

2.5  Summary 

The  use  of  a  reference  image  as  the  a  priori  information  to  guide  the  reduced- 
encoding  data  acquisition  process  does  not  achieve  the  goal  of  exploiting  the  desirable 
truncation  properties  of  the  non- Fourier  basis  sets  and,  at  times,  can  create  dangerous 
artifacts. 
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CHAPTER  3 


IMAGE  RECONSTRUCTION 

The  image  reconstruction  question  is:  given  the  available  data, 

^dyn(^)  —  (-^dynj  ^(^))  ^  ^  •^^yn;  (^•^) 

how  can  the  original  object  function  /dyn  be  recovered?  For  the  reduced-encoding  case, 
the  reconstructed  object  function  is  not  unique,  because  there  is  not  sufficient  data 
available.  In  this  case,  a  priori  information  can  be  used  to  help  select  a  reconstruction 
from  the  set  of  possible  choices.  If  a  priori  information  is  used,  the  key  issue  is 
how  to  effectively  utilize  it  as  a  constraint  during  reconstruction  [71]  to  recover  the 
missing  dynamic  data  ddynl?^),”  i  Three  image  reconstruction  methods  will 

be  discussed  in  this  section  which  use  different  strategies  to  handle  the  unmeasured 
dynamic  data:  zero-padded  reconstruction,  the  data- replacement  technique  and  the 
generalized-series  method. 

3.1  Zero-Padded  Reconstruction 

With  the  zero-padded  reconstruction  method,  the  unmeasured  dynamic  data  are 
simply  set  to  zero.  More  precisely,  from  the  measured  dynamic  data  diya{m),  m  e  A/dyn, 
a  zero-padded  data  set  is  created  as 

f  ddyn(R^')  c  A/dyn 

m)  =  \  (3.2) 

1 0  else, 

where  the  size  of  the  resulting  data  set  depends  on  the  desired  digital  pixel  size.  This 
data  set  will  then  be  processed  with  the  conventional  reconstruction  method. 
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For  example,  with  Fourier-encoded  data,  the  image  can  be  reconstructed  using 
the  discrete  Fourier  transform  as 

M/2-1 

/dy„(nAr)  =  AA:  ddyn(mAA:)e‘^’^^,  (3.3) 

Tn=— M/2 

where  ddyn  is  the  zero-padded  Fourier  encoded  data  set  and  M  is  the  size  of  the  zero- 
padded  data  set.  The  discrete  Fourier  transform  reconstruction  can  be  efficiently 
performed  using  the  Fast  Fourier  transform  (FFT),  which  is  one  of  the  reasons  for 
the  popularity  of  the  Fourier  encoding  method  in  MRI. 

For  orthogonal  wavelet  encoding,  the  original  object  can  be  reconstructed  via  the 
J  level  inverse  wavelet  transform  as 

M2--’ -J  M2-J  _____ 

/(!■)  =  E  2-'<6(2-'r  -  ):)  +  E  E  2V(2'r  -  I),  (3.4) 

fc=l  j=-l  1=1 

where  is  the  zero-padded  set  of  approximation  coefficients  and  is 

the  zero-padded  set  of  wavelet  coefficients.  For  biorthogonal  wavelets,  (j)  and  ip  in 
Eq.  (3.4)  will  be  replaced  by  ^  and  0,  respectively. 

In  the  case  of  SVD  encoding,  the  SVD  synthesis  procedure  is  used  for  image 
reconstruction.  If  the  dynamic  data  were  acquired  using  pure  vertical  SVD  encoding 
as  in  Eq.  (2.11),  the  synthesis  equation  is 

Idyn  =  U  Sdyn,  (3.5) 

where  U  is  the  matrix  containing  the  full  set  of  reference-left  singular  vectors  and 
5dyn  is  the  zero-padded  data  matrix.  To  perform  the  zero-padding  for  this  example, 
M  —  N  rows  of  zeros  are  appended  to  the  bottom  of  the  measured  data  matrix  S'dyn 
of  Eq.  (2.11),  where  M  and  N  are  the  number  of  SVD  encodings  measured  for  the 
reference  and  dynamic  data  sets,  respectively.  The  zero-padding  procedure  can  also 
be  implicitly  handled  using  the  truncated  SVD  reconstruction  procedure.  For  the 
pure  vertical  SVD  encoding  data  set  of  Eq.  (2.11),  the  measured  data  are  multiplied 
by  tZ/sr  as 

fdyn  ~  fZ;V  ‘Fdyni  (3-6) 
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where  Uj^  is  the  matrix  constructed  from  the  A^-left  singular  vectors  selected  for  the 
dynamic  encoding. 


3.2  Data- Replacement  Technique 


With  zero-padded  reconstruction,  it  is  assumed  that  the  unmeasured  data  is  unim¬ 
portant  and  can,  therefore,  be  set  to  zero.  For  the  keyhole  or  data-replacement  tech¬ 
nique  (DRT)  [74, 75],  the  underlying  assumption  is  that  the  dynamic  changes  have 
a  negligible  effect  on  the  unmeasured  dynamic  data  points.  With  this  assumption, 
the  reference  data  can  be  used  to  directly  substitute  for  the  unmeasured  dynamic 
data.  If  the  measured  dynamic  data  are  ddyn(’7i),m  €  and  the  reference  data 

are  dref(m),  m  e  Nre.u  a  full  data  set  is  created  by 


ddyn  (^) 


{ddyn(^)  ^  ^  -^/dyn 
drei{m)  m  G  Wref,  rn  ^  A/dyn- 


(3.7) 


After  the  merged  data  set  is  created,  the  conventional  image  reconstruction  procedure 
is  used  to  reconstruct  the  dynamic  image. 


3.3  Generalized- Series  Technique 

Another  method  to  incorporate  a  priori  reference  information  into  the  reconstruc¬ 
tion  process  is  based  on  the  generalized-series  model.  This  method,  called  the  reduced- 
encoding  imaging  by  generalized-series  reconstruction  (RIGR)  method  [72,73],  is  use¬ 
ful  for  processing  Fourier-encoded  data.  For  wavelet  and  SVD-encoded  data,  there  is 
no  known  better  way  to  use  the  reference  information  in  the  reconstruction  step  than 
the  keyhole  method. 

The  generalized-series  (GS)  model  can  be  described  by 

hy.{r)  =  C{T)  Y.  (3.8) 

^eA/’dyn 

where  C(r)  is  the  constraint  function  and  Afc  is  the  Fourier  encoding  step  in  A;-space. 
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Therefore,  the  RIGR  basis  functions  are  the  set  of  constrained  exponentials 

^{r)  =  C{r)  (3.9) 

which  contain  a  priori  information  from  the  constraint  function.  This  results  in  a  more 
rapidly  convergent  model  than  is  possible  with  the  infinite  complex  exponentials  of 
the  Fourier  series.  The  specific  form  of  the  constraint  function  C{r)  depends  upon 
whether  or  not  a  valid  phase  constraint  is  available  [73].  If  a  valid  phase  constraint 
is  not  available,  the  model  will  be  more  stable  if  the  phase  is  removed,  resulting  in 
basis  functions  of  the  form 

=  (3.10) 

However,  if  a  valid  phase  constraint  is  available,  the  dynamic  changes  will  be  better 

reproduced  if  the  following  basis  functions  are  used 

=  (3.11) 

where  9{r)  is  the  phase  constraint. 

During  the  extrapolation,  data  consistency  between  the  reconstructed  image  and 
the  measured  dynamic  data  is  enforced  during  the  determination  of  the  GS  coeffi¬ 
cients.  If  no  phase  is  used,  the  data  consistency  constraint  will  force  the  GS  parame¬ 
ters  to  absorb  both  the  dynamic  contrast  changes  and  the  phase  variations.  Therefore, 
it  is  desirable  to  have  the  sampling  asymmetry  small  so  that  high-frequency  phase 
variations  will  not  introduce  large  asymmetric  truncation  artifacts.  In  this  case,  the 
N  model  parameters  must  satisfy 
A72-1 

ddyn(?7^)  =  ^  Cn  —  Tl)  —  N/2  <m<  N/2  —  1,  (3.12) 

n=-N/2 

where  symmetric  sampling  of  Arspace  is  assumed  and 

yv  TOO 

-  n)  =  /  |4ef(r)|e-'2-(— (3.13) 

J-oo 

If  phase  constraints  are  used,  the  GS  coefficients  are  forced  to  have  Hermitian  symme¬ 
try  (i.e.,  Cn  =  cl„),  because  the  parameters  need  to  model  only  the  dynamic  contrast 
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changes.  Due  to  this  forced  Hermitian  symmetry,  the  GS  model  can  be  symmetric 
irrespective  of  the  symmetry  of  the  sampled  dynamic  data,  resulting  in  a  higher  order 
model  and,  therefore,  reduced  truncation  artifacts.  The  result  is  two  sets  of  linear 
equations  that  specify  the  GS  parameters  as 

N~No-l 

ddyn  (^)  ~  E  Cndrefim  -  n)  and  (3.14a) 

n=-{N-No-l) 

N-No-l 

ddyni‘>^)  =  -n)  -  No  <m<N  -No  -  1,  (3.14b) 

n=-{N-No-l) 

where  dref  is  the  reference  data.  Because  there  are  more  equations  than  unknowns, 
this  system  should  be  solved  in  the  least  squares  sense.  Substituting  the  coefficients 
of  Eq.  (3.12)  or  (3.14)  into  the  GS  model  will  yield  the  desired  dynamic  image. 


3.4  Generalized  Point  Spread  Function 


As  discussed  in  Chapter  1,  the  width  of  the  point  spread  function  (PSF)  gives 
a  measure  of  the  spatial  resolution  obtained  from  an  imaging  process.  A  common 
definition  of  the  PSF  width  is  the  equivalent  rectangle.  Specifically,  if  h{r)  is  the 
PSF,  the  width  w  of  the  equivalent  rectangle  is  defined  by 

[  /2.(0)rect  f— 'j  dr  =  [  h{r)dr,  (3.15) 

J—oo  \W  /  J—oo 


1^1  <5 


where 

rect(r) 

0  else. 

In  the  case  of  Fourier  series  reconstruction,  the  PSF  is 


-C 


(3.16) 


h(r)  =  Ak 


sm{‘7TNeAkr) 

sin(7rAA:r) 


—inAkr 

^  j 


(3.17) 


where  Ng  is  the  number  of  A:-space  data  points  and  Ak  is  the  step  size  in  krsp&ce,  and 
the  corresponding  PSF  width  Ar  is 


Ar  = 


1 

NgAk' 


(3.18) 
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Note  that  the  width  of  the  PSF  for  Fourier  imaging  depends  only  upon  the  number 
of  encodings  that  are  acquired,  given  a  fixed  Ai-space  step  size.  On  the  other  hand, 
for  the  case  of  constrained  reconstruction  methods  and  specifically  for  the  case  of  GS 
reconstruction,  the  PSF  width  depends  on  both  the  number  of  encodings  and  the  a 
priori  information. 

For  this  reason,  a  generalized  PSF  analysis  was  performed  where  the  reference 
image  was  a  boxcar  function,  and  the  dynamic  change  between  the  reference  and 
dynamic  images  was  a  point  change.  The  location  of  the  point  change  was  varied 
within  the  boxcar.  It  was  either  located  in  the  center,  shifted  one-fourth  of  the  width 
of  the  boxcar  from  the  center,  or  shifted  almost  one-half  (0.49)  of  the  width  of  the 
boxcar  from  the  center.  The  analysis  was  repeated  for  varying  reference  boxcar  widths 
and  a  varying  number  of  dynamic  encodings. 

Example  profiles  are  shown  in  Fig.  3.1  for  the  case  of  eight  dynamic  encodings  and 
a  reference  boxcar  width  of  0.03125  (FOV=l).  To  better  show  details  of  the  plots,  only 
the  center  fourth  of  the  FOV  is  shown  for  each  profile.  Rows  1-3  depict  the  situation 
for  a  centered  point  change  and  a  point  change  shifted  one-fourth  or  one-half  the 
width  of  the  reference  boxcar  from  the  center,  respectively.  Figures  3.1(a)-(c)  are  the 
reference,  dynamic  and  difference  (point  change)  images,  respectively,  reconstructed 
using  512  phase  encodings.  Plots  (d)-(e)  show  the  point  change  reconstructed  using 
the  Fourier  series  and  the  generalized  series,  respectively.  The  reduced  PSF  width  of 
the  generalized  series  reconstruction  when  compared  to  the  Fourier  series  is  due  to  the 
fact  that  the  generalized  series  basis  functions  contain  information  from  the  reference 
image,  so  they  can  more  effectively  reproduce  the  dynamic  image  for  a  given  number 
of  encodings  than  the  infinite  complex  exponentials  of  the  Fourier  series.  Note  that 
in  this  case,  the  reference  image  contains  no  edge  information  for  the  new  dynamic 
feature. 

The  effects  of  the  reference  boxcar  width  and  number  of  dynamic  encodings  on 
the  width  of  the  PSF  can  be  seen  in  Fig.  3.2.  As  would  be  expected,  the  width  of  the 
GS  PSF  decreases  with  an  increasing  number  of  dynamic  encodings.  This  also  holds 
true  for  the  Fourier  series  reconstruction.  Note  also  that  the  width  of  the  GS  PSF 
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(1) 


(2) 


(3) 


Figure  3.1:  PSF  Profiles  I:  Rows  1-3  show  the  PSF  results  for  a  delta 
function  change  that  is  centered  in  the  reference  boxcar,  shifted  by  one- 
fourth  the  width  of  the  reference  boxcar,  and  shift  by  just  under  one-half 
(0.49)  of  the  width  of  the  reference  boxcar,  respectively.  The  width  of 
the  reference  boxcar  was  0.03125  (FOV=l),  but  only  the  center  fourth  of 
the  plot  is  shown  for  better  visualization,  (a)-(c)  The  reference  image, 
dynamic  image  and  point  change  image  (difference  image),  respectively, 
reconstructed  using  512  phase  encodings,  (d)-(e)  The  point  change  recon¬ 
structed  using  the  Fourier  series  and  generalized  series,  respectively,  using 
eight  dynamic  encodings.  Note  that  (a)-(c)  are  on  a  different  scale  than 
(d)-(e). 
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PSF  vs.  Number  of  Dynamic  Encodings 
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Figure  3.2:  PSF  vs  Number  of  Dynamic  Encodings:  This  plot  shows  the 
relationship  of  the  width  of  the  PSF  versus  the  number  of  dynamic  encod¬ 
ings  used  for  the  Fourier  series  and  generalized  series.  In  the  simulations 
used  to  generate  this  plot,  the  point  change  was  centered  in  the  refer¬ 
ence  boxcar.  Note  the  reduced  PSF  width  of  the  GS  as  compared  to  the 
Fourier  series  for  all  reference  boxcar  widths  and  all  numbers  of  dynamic 
encodings. 


Centered  Point  Change 


decreases  for  decreasing  width  of  the  reference  boxcar.  Although  no  additional  direct 
information  is  available  about  the  dynamic  change,  the  reference  information  of  the 
narrower  boxcars  is  better  able  to  constrain  the  reconstruction  of  the  new  dynamic 
feature.  For  every  reference  boxcar  width,  the  width  of  the  PSF  is  smaller  for  the 
generalized  series  than  for  the  Fourier  series.  In  the  limit  that  the  reference  boxcar 
width  equals  the  FOV,  the  generalized  series  has  the  same  PSF  width  as  the  Fourier 
series.  In  fact,  for  this  case  of  no  effective  a  priori  knowledge,  the  generalized  series 
reduces  to  the  Fourier  series.  This  is  desirable,  because  in  the  case  of  no  a  priori 
knowledge,  the  Fourier  series  would  be  the  optimal  reconstruction  method  (in  the 
least  squares  sense). 

The  GS  PSF  width  depends  not  only  on  the  reference  boxcar  width,  but  also  on 
the  proximity  of  the  point  change  to  the  nearest  edge.  This  is  illustrated  in  Fig.  3.3, 
which  shows  the  PSF  width  as  a  function  of  the  width  of  the  reference  boxcar  for  the 
Fourier  series  and  the  generalized  series  with  various  locations  of  the  point  change. 
Note  that  the  PSF  width  of  the  generalized  series  decreases  with  increased  proximity 
to  an  edge  of  the  reference  boxcar  (i.e.,  increasing  shift  from  the  center).  This  occurs 
because  the  closer  edge  enables  the  GS  basis  functions  to  more  effectively  constrain 
the  dynamic  change  reconstruction. 

If  the  high-resolution  reference  image  contains  edge  information  for  the  dynamic 
feature,  the  GS  PSF  will  be  further  improved  as  illustrated  in  Fig.  3.4.  As  before, 
rows  1-3  show  the  reconstruction  of  a  dynamic  point  change  that  is  centered  in  the 
reference  boxcar  and  shifted  one-fourth  or  approximately  one-half  (0.49)  the  width 
of  the  reference  boxcar  from  the  center,  respectively.  Figures  3.4 (a)- (d)  show  the 
baseline  reference  image,  the  active  reference  image,  the  dynamic  image,  and  the 
point  change  image  (difference  between  the  dynamic  image  and  the  baseline  reference 
image),  respectively,  reconstructed  using  512  phase  encodings.  Plot  (e)  shows  the 
reconstruction  of  the  point  change  obtained  using  the  GS  model  with  the  active 
reference  image  and  eight  dynamic  encodings.  It  is  easy  to  see  the  improved  PSF  by 
comparison  to  the  corresponding  rows  of  Fig.  3.1(e). 
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PSF  vs.  Width  of  Reference  Boxcar 
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Figure  3.3:  PSF  vs  Width  of  Reference  Boxcar:  This  plot  shows  the  re¬ 
lationship  between  the  width  of  the  PSF  and  the  width  of  the  reference 
boxcar  for  the  Fourier  series  and  generalized  series  for  three  locations  of 
the  point  change.  The  simulations  to  generate  this  plot  used  64  dynamic 
encodings. 


32 


(1) 


(2) 


(3) 


Figure  3.4:  PSF  Profiles  II:  Rows  1-3  show  the  PSF  results  for  a  delta 
function  change  that  is  centered  in  the  reference  boxcar,  shifted  by  one- 
fourth  the  width  of  the  reference  boxcar,  and  shift  by  just  under  one-half 
(0.49)  of  the  width  of  the  reference  boxcar,  respectively.  The  width  of  the 
reference  boxcar  was  0.03125  (FOV=l),  but  only  the  center  fourth  of  the 
plot  is  shown  for  better  visualization,  (a)-(d)  The  baseline  reference  im¬ 
age,  the  active  reference  image,  the  dynamic  image,  and  the  point  change 
image,  respectively,  reconstructed  using  512  phase  encodings,  (e)  The 
point  change  reconstructed  using  the  GS  model  with  the  active  reference 
image  and  eight  dynjunic  encodings.  Note  that  (a)-(c)  are  on  a  different 
scale  than  (d)-(e). 
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In  summary,  the  GS  model  has  a  smaller  PSF  width  than  the  Fourier  series  even 
if  the  reference  image  contains  no  edge  information  for  the  dynamic  feature  because 
nearby  edges  help  constrain  the  reconstruction  of  the  dynamic  feature.  If  the  reference 
image  does  contain  edge  information  about  the  dynamic  feature,  the  PSF  width  will 
be  further  reduced. 


3.5  Discussion 


For  image  reconstruction,  the  use  of  a  priori  information  in  the  constrained  re¬ 
construction  techniques  such  as  keyhole  and  RIGR  can  improve  the  resulting  image 
appearance  over  that  obtained  through  simple  zero-padded  reconstruction.  However, 
the  actual  improvement  gained  will  vary  greatly  depending  upon  the  method  of  uti¬ 
lizing  the  a  priori  information.  Although  the  keyhole  method  has  been  applied  to 
reduced-encoding  dynamic  imaging  applications  due  to  its  simplicity  [67,81,82],  it 
is  important  to  note  that  the  Fourier-keyhole  method  can  only  track  the  dynamic 
changes  at  the  low  resolution  of  1/A^AA:,  where  N  is  the  number  of  reduced  dynamic 
encodings^  [81,83,84].  This  occurs  because  the  keyhole  method  has  the  same  PSF 
as  the  zero-padded  Fourier  series  for  reconstructing  the  dynamic  changes.  This  can 
easily  be  seen  as 


/diff(r) 


•fdyn(^)  -frefC^) 

^  Y.  ddyn(n)e-'2-^"^  ^ 

X)  -  dref{n)  }  . 


(3.19a) 

(3.19b) 

(3.19c) 

(3.19d) 


On  the  other  hand,  the  spatial  resolution  of  the  RIGR  image  depends  upon  both  the 
dynamic  and  reference  data  sets.  If  well-defined  boundaries  for  the  feature  exist  in  the 
reference  image,  it  will  be  reconstructed  with  the  spatial  resolution  of  the  reference 

^For  the  wavelet-keyhole  and  SVD-keyhole  methods,  the  resolution  at  which  the  dynamic  changes 
can  be  followed  will  depend  on  the  relationship  between  the  dynamic  changes  and  the  selected  wavelet 
or  SVD  encoding  profiles. 
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image  1/MA/c,  where  M  is  the  number  of  reference  encodings  [73].  However,  as 
was  shown  in  the  previous  section,  even  if  the  dynamic  changes  introduce  new  edges 
to  the  image,  the  spatial  resolution  of  these  features  can  be  improved  over  1/A^AA:, 
because  nearby  edges  in  the  reference  image  help  to  constrain  the  reconstruction  of 
the  new  features.  Therefore,  with  the  same  number  of  dynamic  encodings,  RIGR  can 
reconstruct  the  dynamic  changes  with  a  greater  spatial  resolution  than  is  possible  with 
Fourier-keyhole.  The  effect  of  this  can  be  seen  in  Fig.  3.5  in  which  (a)-(b)  are  the 
reference  and  dynamic  images,  respectively,  reconstructed  using  256  phase  encodings. 
Figure  3.5(c)  is  the  difference  between  these  two  images,  effectively  an  image  of  the 
dynamic  changes.  Figures  3.5(d)-(e)  show  the  difference  image  reconstructed  with 
32  dynamic  encodings  using  Fourier-keyhole  (or,  equivalently,  zero-padded  Fourier 
series)  and  RIGR,  respectively.  Note  that  the  RIGR  method  follows  the  dynamic 
changes  at  a  much  higher  resolution  than  Fourier-keyhole. 

In  addition,  image  artifacts  can  arise  due  to  data  inconsistency  between  the  refer¬ 
ence  and  dynamic  data  sets  with  the  keyhole  method  [85,86].  In  the  keyhole  method 
where  the  reference  data  are  simply  pasted  on  the  dynamic  data,  there  is  no  guar¬ 
antee  of  consistency  between  the  reference  and  dynamic  data  sets.  This  can  lead  to 
image  artifacts.  In  contrast,  the  unmeasured  high-frequency  dynamic  data  extrapo¬ 
lated  using  the  GS  model  are  (TV  —  1)‘^  order  continuous  with  the  measured  dynamic 
data  [73].  This  leads  to  reduced  artifacts  in  the  reconstructed  dynamic  image.  This 
is  illustrated  in  Fig.  3.6  in  w^hich  (a)-(b)  are  the  reference  and  dynamic  images,  re¬ 
spectively,  reconstructed  with  128  phase  encodings.  Images  (c)-(d)  show  the  dynamic 
image  reconstructed  with  16  dynamic  encodings  using  the  Fourier-keyhole  and  RIGR 
methods,  respectively.  Note  the  data  inconsistency  artifacts  in  the  Fourier-keyhole 
image,  which  are  equivalent  to  a  high-pass  filtered  version  of  the  reference  image.^ 
On  the  other  hand,  the  RIGR  method  can  handle  the  data  inconsistency  between  the 
reference  and  dynamic  data  sets  due  to  the  GS  model-based  extrapolation. 

^Data  inconsistency  artifacts  will  also  arise  with  the  wavelet-keyhole  and  SVD-keyhole  reconstruc¬ 
tion  methods,  although  the  appearance  of  the  artifacts  will  differ  from  those  of  the  Fourier-keyhole 
method. 
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Figure  3.5:  Dynamic  Change  Images:  (a)-(c)  The  reference  and  dynamic 
images  and  the  difference  between  the  two,  respectively,  (d)-(e)  The  dif¬ 
ference  image  reconstructed  using  32  dynamic  encodings  with  Fourier- 
keyhole  (or,  equivalently,  zero-padded  Fourier  series)  and  RIGR,  respec¬ 
tively.  Note  the  improved  reconstruction  of  the  dynamic  changes  with  the 
RIGR  method. 
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Figure  3.6:  Keyhole  Data  Inconsistency  Artifact:  (a)-(b)  The  reference 
and  dynamic  images,  respectively,  reconstructed  using  128  phase  encod¬ 
ings.  (c)-(d)  The  dynamic  image  reconstructed  with  16  dynamic  encod¬ 
ings  using  Fourier-keyhole  and  RIGR,  respectively.  The  phase  encoding 
direction  is  horizontal.  Note  the  edge  artifacts  that  appear  in  the  Fourier- 
keyhole  image. 


36 


3.6  Summary 


If  the  a  priori  information  that  is  available  is  a  reference  image,  it  is  better  used 
in  the  image  reconstruction  process  than  the  data-acquisition  step.  For  image  recon¬ 
struction  with  Fourier  encoding,  the  RIGR  method  is  the  best  way  to  extrapolate  the 
unmeasured  data  using  the  a  priori  constraints  due  to  the  higher  resolution  tracking 
of  the  dynamic  changes  and  the  reduced  data  inconsistency  artifacts. 
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CHAPTER  4 


TWO-REFERENCE  RIGR 


With  the  Fourier  series,  the  spatial  resolution  of  the  reconstructed  image  depends 
solely  upon  the  number  of  terms  used  in  the  series.  On  the  other  hand,  the  spa¬ 
tial  resolution  obtained  with  the  generalized  series  (GS)  model  depends  on  both  the 
number  of  terms  and  the  GS  basis  functions.  Due  to  the  limited  nQmber  of  terms 
in  the  GS  model,  it  would  be  better  if  only  the  areas  of  change  were  represented  in 
the  GS  basis  functions.  The  proposed  two-reference  reduced-encoding  imaging  by 
generalized-series  reconstruction  (TRIGR)  method  [87,88]  tries  to  achieve  this  by 
suppressing  the  background  information  in  the  GS  basis  functions  through  the  use  of 
a  second  high-resolution  active  reference  image. 

The  motivation  behind  the  research  presented  here  is  the  consideration  that,  in 
some  dynamic  imaging  applications,  the  dynamic  process  may  change  the  appearance 
of  parts  of  the  image,  while  other  parts  remain  relatively  unchanged.  For  example, 
in  contrast-enhanced  dynamic  imaging  of  breast  cancer,  where  the  aim  is  to  track 
the  changes  that  occur  in  the  breast  for  several  minutes  following  the  injection  of  a 
contrast  agent,  the  contrast  between  the  tumor  and  the  normal  tissue  may  change 
drastically  over  the  dynamic  imaging  period.  The  TRIGR  method  is  well-suited  for 
this  application,  because  a  second  high-resolution  reference  image  can  be  acquired 
following  the  dynamic  imaging  period  when  the  contrast  agent  is  strongly  visible 
in  the  image.  The  second  reference  image  can  be  used  to  suppress  the  background 
information  in  the  generalized-series  basis  functions,  thus  improving  the  reconstructed 
dynamic  image. 
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4.1  Two- Reference  RIGR 

Data  acquisition  for  the  proposed  method  is  similar  to  the  original  RIGR  method. 
A  high-resolution  baseline  reference  data  set  is  acquired  where  the  number  of  phase 
encodings  is  chosen  to  satisfy  the  spatial  resolution  requirements.  This  is  followed  by 
the  acquisition  of  a  series  of  reduced-encoding  dynamic  data  sets  where  the  number  of 
Fourier  encodings  per  set  is  chosen  to  give  the  desired  temporal  resolution.  The  high- 
resolution  active  reference  data  set  is  typically  acquired  at  the  end  of  the  dynamic 
imaging  period,  although  for  some  applications  it  may  be  desirable  to  place  it  at 
another  point  in  the  experimental  procedure. 

The  proposed  method  suppresses  background  information  in  the  GS  basis  func¬ 
tions  through  the  use  of  a  difference  reference  image  which  is  created  as 

/ref  ~  •^■{dref}  —  •^'{dactive  ^baseline}) 

where  ^active  and  dbaseiine  are  the  high-resolution  active  and  baseline  reference  data 
sets,  respectively.  The  GS  basis  functions  become  the  set  of  constrained  exponentials 
resulting  in  the  model 

N/2-1 

/dil(r)  =  4f(r)  i:  (4.2) 

n=-N/2 

where  N  is  the  number  of  dynamic  encodings  and  symmetric  A:-space  sampling  has 
been  assumed.  Note  that  the  TRIGR  method  directly  reconstructs  an  image  of  the 
dynamic  change.  As  such,  the  difference  between  the  reduced  dynamic  encodings  and 
the  corresponding  baseline  reference  encodings  are  used  to  determine  the  GS  model 
coefficients.  Specifically,  the  GS  model  coefficients  c„  are  obtained  through  a  fitting 
step  to  maintain  data  consistency  as 

N/2-1 

ddArn)=  Cndrefim-n),  (4.3) 

n=-JV/2 

where  the  dynamic  difference  data  are  obtained  by  subtracting  from  the  dynamic 
data  the  corresponding  encodings  of  the  baseline  reference  image  as 

-N  N 

ddisij^^k')  —  ddyniP'^k')  0?baseline(^^^)  ^  ^  !•  (‘^•^) 
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Plugging  these  coefficients  into  Eq.  (4.2)  will  yield  the  reconstructed  dynamic  change 
image.  If  the  dynamic  image  itself  is  desired,  it  can  be  generated  by  adding  the 
complex  dynamic  change  image  to  the  baseline  reference  image  as 

-fdyii(^)  ~  .fbaseline(^)  4”  -fdifF(^))  (^■^) 

where  /baseline(^)  is  reconstructed  from  dbaseiine(?^^^)  using  the  standard  Fourier  tech¬ 
nique. 

The  background  information  suppression  achieved  with  this  technique  can  result 
in  significant  improvement  in  the  generalized  PSF,  as  illustrated  in  Fig.  4.1.  As  be¬ 
fore,  rows  1-3  show  the  reconstruction  of  a  dynamic  point  change  that  is  centered 
in  the  reference  boxcar  and  shifted  one-fourth  or  approximately  one-half  (0.49)  the 
width  of  the  reference  boxcar  from  the  center,  respectively.  Figures  4.1(a)-(d)  show 
the  baseline  reference  image,  the  active  reference  image,  the  dynamic  image,  and  the 
point  change  image  (difference  between  the  dynamic  image  and  the  baseline  refer¬ 
ence  image),  respectively,  reconstructed  using  512  phase  encodings.  Plot  (e)  shows 
the  reconstruction  obtained  using  the  GS  model  with  the  difference  reference  image 
(active  -  baseline)  and  eight  dynamic  encodings.  It  is  easy  to  see  the  improvement 
in  the  PSF  due  to  the  suppression  of  the  background  information  in  the  difference 
reference  image  by  comparison  to  the  corresponding  rows  of  Fig.  3.4(e).  Although 
the  case  shown  in  Fig.  3.4  has  the  same  amount  of  edge  information  for  the  dynamic 
feature,  the  background  suppression  provides  additional  improvement.  Of  course,  this 
is  the  ideal  case  of  complete  background  suppression.  As  the  background  information 
is  less  completely  suppressed,  the  improvement  seen  with  the  TRIGR  method  will 
not  be  as  great. 

One  could  consider  employing  a  similar  methodology  with  the  keyhole  constrained 
reconstruction  method  by  appending  the  high-frequency  difference  reference  data  to 
the  low-frequency  dynamic  difference  data  sets  followed  by  the  inverse  Fourier  trans¬ 
form.  Although  the  resulting  dynamic  images  will  be  high-resolution,  it  is  easy  to 
prove  that  the  actual  dynamic  signal  changes  will  still  be  reconstructed  with  low- 
resolution.  This  behavior  is  similar  to  the  single- reference  image  case  analyzed  by 
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(1) 


(2) 


(3) 


Figure  4.1:  PSF  Profiles  III:  Rows  1-3  show  the  PSF  results  for  a  delta 
function  change  that  is  centered  in  the  reference  boxcar,  shifted  by  one- 
fourth  the  width  of  the  reference  boxcar,  and  shift  by  just  under  one-half 
(0.49)  of  the  width  of  the  reference  boxcar,  respectively.  The  width  of 
the  reference  boxcar  was  0.03125  (FOV=l),  but  only  the  center  fourth 
of  the  plot  is  shown  for  better  visualization,  (a)-(d)  The  baseline  refer¬ 
ence  image,  the  active  reference  image,  the  dynamic  image,  and  the  point 
change  image,  respectively,  reconstructed  using  512  phase  encodings,  (e) 
The  point  change  reconstructed  using  the  GS  model  with  the  a  difference 
reference  image  (active-baseline)  and  eight  dynamic  encodings.  Note  that 
(a)-(c)  are  on  a  different  scale  than  (d)-(e). 
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Spraggins  [84]  and  Hu  [83];  that  is,  no  benefit  is  gained  from  the  use  of  two-reference 
images  in  this  keyhole  scheme. 

4.2  Results  and  Discussion 

The  proposed  method  was  compared  to  the  zero-padded  Fourier  transform  and 
the  original  RIGR  method  using  computer  simulations.  From  a  sequence  of  high- 
resolution  data  sets,  baseline  and  active  reference  data  sets  were  chosen.  The  central 
N  encodings  of  the  remaining  data  sets  were  used  as  the  dynamic  data  sets,  and 
the  dynamic  image  was  reconstructed  with  the  three  methods.  The  high-resolution 
dynamic  image  was  reconstructed  for  comparison. 

As  mentioned  before,  the  improvement  of  the  proposed  technique  over  the  original 
RIGR  method  depends  upon  the  degree  of  background  suppression.  The  ideal  case 
is  if  the  difference  image  completely  suppresses  the  background  information.  This 
case  is  illustrated  in  Figs.  4.2  rows  1-3  which  show  the  resulting  images  and  profiles 
through  the  upper  and  lower  set  of  lesions,  respectively,  for  a  contrast-enhanced 
simulation  in  which  the  background  tissue  is  completely  suppressed  in  the  difference 
reference  image.  In  these  figures,  (a)-(c)  are  the  precontrast  reference,  postcontrast 
reference  and  dynamic  images,  respectively,  reconstructed  using  128  phase  encodings. 
Images  (d)-(e)  were  reconstructed  using  RIGR  with  16  dynamic  encodings  and  the 
precontrast  or  postcontrast  reference  image,  respectively.  Image  (f)  was  reconstructed 
using  16  dynamic  encodings  with  the  proposed  method.  The  TRIGR  method  results 
in  an  improved  reconstruction  of  the  dynamic  image  because  the  GS  basis  functions 
need  only  represent  the  dynamic  changes. 

As  the  difference  reference  image  suppresses  the  background  less  effectively,  the 
performance  improvement  with  the  proposed  method  decreases.  This  is  demonstrated 
in  Fig.  4.3  rows  1-3  which  show  the  images  and  profiles  through  the  upper  and  lower 
set  of  lesions,  respectively,  resulting  from  a  simulation  in  which  the  background  tissue 
is  not  suppressed  completely  in  the  difference  reference  image.  Figures  4.3  (a)- (c)  are 
the  precontrast  reference,  postcontrast  reference  and  dynamic  images,  respectively. 
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Figure  4.2:  TRIGR  With  Complete  Background  Suppression:  Images  (row 
1)  and  profiles  through  the  upper  and  lower  set  of  lesions  (rows  2-3,  re¬ 
spectively).  (a)-(c)  The  baseline  reference,  active  reference  and  dynamic 
images,  respectively,  reconstructed  using  128  phase  encodings,  (d)-(e)  The 
dynamic  image  reconstructed  using  16  dynamic  encodings  with  RIGR  us¬ 
ing  the  baseline  and  active  reference  images,  respectively,  (f)  The  dynamic 
image  reconstructed  using  16  dynamic  encodings  with  TRIGR. 

reconstructed  using  128  phase  encodings.  Images  (c)-(f)  were  reconstructed  with 
16  dynamic  encodings  using  RIGR  with  the  precontrast  reference,  RIGR  with  the 
postcontrast  reference  and  the  proposed  method,  respectively.  The  improvement 
from  using  TRIGR  is  much  less  in  this  case  than  in  Fig.  4.2,  because  the  background 
is  not  as  well  suppressed. 

The  results  of  all  constrained  image  reconstruction  techniques  depend  upon  the 
validity  of  the  a  priori  constraints.  As  with  RIGR,  the  results  of  the  proposed  method 
depend  upon  the  choice  of  the  baseline  reference  image.  In  addition,  the  active 
reference  image  will  also  affect  the  outcome  of  the  TRIGR  reconstruction  process. 
As  would  be  expected,  the  more  similar  the  baseline  or  active  reference  image  is  to 
the  dynamic  image,  the  better  the  reconstruction.  This  effect  is  illustrated  in  the 
following  simulation,  which  investigates  the  performance  of  the  TRIGR  algorithm 
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Figure  4.3:  TRIGR  With  Incomplete  Background  Suppression:  Images 
(row  1)  and  profiles  through  the  upper  and  lower  set  of  lesions  (rows 
2-3,  respectively),  (a)-(c)  The  baseline  reference,  active  reference  and 
dynamic  images,  respectively,  reconstructed  using  128  phase  encodings. 
(d)-(e)  The  dynamic  image  reconstructed  using  16  dynamic  encodings  with 
RIGR  using  the  baseline  and  active  reference  images,  respectively,  (f)  The 
dynamic  image  reconstructed  using  16  dynamic  encodings  with  TRIGR. 

with  different  active  reference  images.  The  reference  images  shown  in  the  top  row 
of  Fig.  4.4  are  from  different  time  points  in  a  contrast-enhanced  dynamic  imaging 
simulation.  Rows  2  and  3  show  profiles  through  the  upper  and  lower  set  of  lesions, 
respectively.  All  of  the  images  in  this  figure  were  reconstructed  using  128  phase 
encodings. 

Figure  4.5  rows  1-3  show  the  resulting  images  and  profiles  through  the  upper  and 
lower  set  of  lesions,  respectively,  of  a  TRIGR  simulation  in  which  Fig.  4.4(d)  was 
used  as  the  dynamic  image.  For  comparison  purposes.  Fig.  4.5(a)  is  the  dynamic 
image  reconstructed  using  128  phase  encodings.  Figures  4.5(b)-(f)  show  the  TRIGR 
reconstruction  obtained  with  eight  dynamic  encodings  using  Fig.  4.4(a)  as  the  base¬ 
line  reference  image  and  Figs.  4.4(b)- (f),  respectively,  as  the  active  reference  image. 
Therefore,  (b)-(c)  use  an  active  reference  that  is  less  enhanced  than  the  dynamic 
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Figure  4.4:  Contrast-enhanced  Simulation  Reference  Images:  Row  1  shows 
images  from  different  time  points  in  a  contrast-enhanced  dynamic  imaging 
simulation  reconstructed  using  128  phase  encodings.  Rows  2  and  3  show 
profiles  through  the  upper  and  lower  set  of  lesions,  respectively. 

image,  (e)-(f)  use  an  active  reference  that  is  more  enhanced  than  the  dynamic  image, 
and  (d)  uses  the  ideal  active  reference  image,  i.e.,  the  dynamic  image  itself.  Note 
that  the  lesion  reconstruction  is  improved  with  a  more  similar  reference  image.  This 
suggests  that,  in  some  situations,  it  may  be  desirable  to  acquire  reference  images 
at  various  points  in  the  experimental  procedure  and  then  use  the  appropriate  two 
references  for  the  reconstruction  of  a  particular  dynamic  image. 

The  effect  of  the  choice  of  active  reference  image  can  be  seen  in  actual  MRI  images 
such  as  those  shown  in  Fig.  4.6.  The  data  for  this  simulation  is  from  a  contrast- 
enhanced  dynamic  imaging  study  of  a  rat  with  breast  cancer.  A  spin-echo  sequence 
(TR  300/TE  20)  was  used  to  collect  the  data.  A  high-resolution  precontrast  reference 
data  set  was  first  acquired.  The  contrast  agent  was  then  injected  and  a  series  of 
dynamic  data  sets  was  acquired  as  the  contrast  agent  washed  into  the  slice.  One  of  the 
data  sets  from  the  dynamic  imaging  period  was  selected  as  the  dynamic  data  set  for 
this  simulation.  Two  of  the  other  data  sets  were  used  for  the  active  reference  images. 
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Figure  4.5:  TRIGR  with  Different  Active  Reference  Images:  Images  (row 
1)  and  profiles  through  the  upper  and  lower  set  of  lesions  (rows  2-3,  re¬ 
spectively).  (a)  Dynamic  image  reconstructed  using  128  phase  encod¬ 
ings.  The  remaining  images  (b)-(f)  were  reconstructed  with  TRIGR  using 
Fig.  4.4(a)  as  the  baseline  reference,  Fig.  4.4(d)  as  the  dynamic  image  and 
Fig.  4.4(b)-(f),  respectively,  as  the  active  reference.  The  TRIGR  images 
were  reconstructed  using  eight  dynamic  phase  encodings. 

Figures  4.6  (a)-(d)  are  the  precontrast  reference  image,  a  less-enhanced  postcontrast 
reference  image,  a  more-enhanced  postcontrast  reference  image,  and  the  dynamic 
image,  respectively,  reconstructed  using  256  phase  encodings.  Figures  4.6(e)-(f)  show 
the  TRIGR  reconstruction  resulting  from  using  eight  dynamic  encodings  with  (b) 
and  (c),  respectively,  as  the  postcontrast  reference  image.  The  influence  of  the  active 
reference  image  can  easily  be  seen  in  (e)-(f).  Because  only  eight  dynamic  encodings 
were  used,  the  effect  of  the  active  reference  image  on  the  reconstructed  dynamic  image 
is  stronger  than  if  a  larger  number  of  dynamic  encodings  were  used. 
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Figure  4.6;  Contrast-enhanced  Study  Images:  (a)-(c)  are  the  precontrast 
reference  image  and  two  different  postcontrast  reference  images  recon¬ 
structed  using  256  phase  encodings,  (d)  The  dynamic  image  reconstructed 
using  256  phase  encodings,  (e)-(f)  The  dynamic  image  reconstructed  us¬ 
ing  TRIGR  with  eight  dynamic  encodings  using  (b)-(c),  respectively,  as 
the  postcontrast  reference  image.  This  corresponds  to  a  less-enhanced 
postcontrast  reference  image  for  (e)  and  a  more-enhanced  postcontrast 
reference  image  for  (f). 

4.3  Summary 

The  proposed  dynamic  imaging  method  can  produce  improved  dynamic  images 
over  the  original  RIGR  method,  if  the  difference  reference  image  can  provide  a  level 
of  suppression  of  the  background  information.  The  benefit  derived  from  the  proposed 
method  increases  with  the  level  of  background  suppression.  The  method  should  prove 
valuable  in  applications  such  as  contrast-enhanced  dynamic  studies,  functional  imag¬ 
ing,  and  interventional  MRI. 
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CHAPTER  5 


GS  WITH  EXPLICIT  BOUNDARY 
CONSTRAINTS 


Recall  that  the  GS  reconstruction  formula  is  given  by 


(5,1) 


Clearly,  the  resulting  image  quality  will  improve  as  the  basis  functions  approach  those 
that  would  be  derived  from  the  dynamic  image  itself,  i.e.,  if 


C{r)e 


-i27TnAkr 


(5.2) 


In  particular,  because  the  GS  coefficients  Cn  are  global,  the  better  the  contrast  re¬ 
lationship  in  the  constraint  function  C(r)  matches  that  in  the  dynamic  image,  the 
better  the  resulting  reconstruction.  If  the  contrast  in  the  reference  and  dynamic  im¬ 
ages  is  the  same,  it  is  helpful  to  build  the  reference  contrast  information  into  the 
GS  basis  functions.  On  the  other  hand,  if  the  contrast  relationship  is  not  the  same, 
only  the  boundary  information  from  the  reference  image  should  be  used  to  constrain 
the  reconstruction.  The  proposed  method  uses  explicit  edge  information  from  the 
reference  image  along  with  contrast  information  from  the  dynamic  data  to  improve 
the  GS  basis  functions.  This  differs  from  the  GS  methods  described  earlier,  which 
use  the  dynamic  data  only  to  determine  the  GS  coefficients. 

The  proposed  method  is  particularly  applicable  for  dynamic  imaging  applications 
such  as  contrast-enhanced  studies  in  which  the  contrast  relationship  between  regions 
of  the  image  can  differ  significantly  from  that  seen  in  a  precontrast  or  postcontrast 
reference  image.  By  modifying  the  contrast  in  the  reference  image  to  more  closely 
align  with  that  of  the  dynamic  image,  the  GS  basis  functions  derived  from  the  new 
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reference  image  will  be  better  able  to  represent  the  dynamic  image,  leading  to  an 
improved  dynamic  image  reconstruction. 


5.1  Generalized  Series  Model  With  Explicit  Edge 
Constraints 


The  GS  model  represents  the  dynamic  image  as  a  contrast  modulation  of  the 
underlying  high-resolution  reference  image  as 


/dy,{r)  = /„,(r)  Z  c{n)e-“'"“'-  =  /„,(r)C„(r),  (5.3) 

n6A/’(jyn 

where  Cmii")  is  the  contrast  modulation  function.  Therefore,  the  dynamic  changes 
reproduced  by  the  GS  model  may  not  have  the  same  resolution  as  the  reference  image, 
because  the  contrast  modulation  function  is  limited  by  the  number  of  terms  in  the 
model,  which  is  in  turn  determined  by  the  number  of  dynamic  encodings.  To  alleviate 
this  problem,  the  proposed  method  seeks  to  incorporate  dynamic  information  into  the 
basis  functions  in  addition  to  the  GS  parameters. 

To  illustrate  this  idea,  consider  segmenting  the  reference  image  and  selecting 
Nreg  <  iV  of  the  most  important  regions  Tin,  where  N  is  the  number  of  dynamic 
encodings.  To  a  first  approximation,  the  difference  between  the  dynamic  and  refer¬ 
ence  images  can  be  represented  as 


Idisi'f’) 


dji  r  €  'R'n 
0  else. 


(5.4) 


In  the  ID  case,  this  can  also  be  expressed  as 

Nreg 

Im{r)  [U{r  -  e„)  -  U{r  -  e„+i)],  (5.5) 

71=1 

where  U  is  the  unit  step  function  and  Cn  are  the  edges  extracted  from  the  high- 
resolution  reference  image.  The  dn  can  be  obtained  by  fitting  the  difference  between 
the  dynamic  data  and  the  corresponding  encodings  of  the  reference  data  to  the  model 
in  Eq.  (5.5).  The  resulting  dynamic  change  image  will  have  high-resolutioii  dynamic 
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variations,  because  the  region  boundary  locations  are  not  limited  by  the  number  of 
dynamic  encodings.  To  utilize  this  information  in  the  GS  model,  a  new  reference 
image  is  created  as 

/ref(0=4ef(r’)+W0-  (5.6) 

The  resulting  GS  basis  functions  will  contain  dynamic  information. 

Higher-order  dynamic  changes  can  be  modeled  through  the  use  of  a  higher-order 
polynomial  model  or  another  set  of  basis  functions,  such  as  wavelets,  in  each  region. 
However,  this  leads  to  problems  for  this  application  due  to  the  large  amount  of  data 
that  would  be  required  for  reliable  fitting.  For  the  example  above  of  modeling  each 
region  with  a  boxcar  function,  if  the  edges  are  known  exactly,  one  parameter  per 
region  is  necessary  to  characterize  the  function  exactly.  Higher-order  behavior  would 
require  a  corresponding  increase  in  the  number  of  parameters.  Due  to  the  small 
amount  of  dynamic  data  that  is  available  in  this  application,  the  number  of  regions 
and  the  model  order  in  each  region  must  be  limited  to  ensure  the  stability  of  the 
fitting  step.  This  leads  to  an  averaging  of  adjacent  regions,  resulting  in  a  blurring 
effect. 

To  alleviate  this  problem,  the  proposed  method  determines  the  average  signal 
magnitude  change  for  each  of  the  regions  detected  by  the  edge  extraction.  Therefore, 
because  no  fitting  step  is  required,  all  of  the  regions  can  be  used,  which  reduces  the 
averaging  effect.  As  mentioned  earlier,  the  extracted  edges  are  used  to  determine 
the  dynamic  change  between  the  reference  image  and  dynamic  image,  as  opposed  to 
the  new  reference  image  itself.  This  will  further  reduce  the  averaging  effect,  as  well 
as  reduce  the  adverse  effects  of  edges  that  are  not  accurately  detected  by  the  edge 
extraction  step. 

For  this  application,  a  multiscale  edge  detection  approach  is  desirable,  because  it 
is  expected  that  image  edges  will  appear  at  many  scales.  In  addition,  the  boundary 
extraction  step  should  result  in  closed  edges  to  allow  easy  determination  of  the  regions 
in  the  image.  The  particular  method  used  for  this  research  was  the  multiresolution 
edge  detection  scheme  developed  by  Ahuja  [89,90].  In  contrast  to  the  multiscale  edge 
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detection  that  can  be  performed  by  some  wavelet  transforms,  the  concept  of  scale 
in  the  multiresolution  edge  detection  scheme  relates  to  both  physical  proximity  and 
greyscale  similarity,  as  opposed  to  the  size  of  the  edge  itself,  which  is  beneficial  for 
this  application. 

The  specific  steps  for  applying  the  proposed  method  to  a  high-resolution  refer¬ 
ence  image  and  a  low-resolution  dynamic  data  set  are  as  follows.  First,  a  zero-padded 
difference  image  /dis  is  created  by  subtracting  the  corresponding  reference  encodings 
from  the  dynamic  encodings  and  reconstructing  using  the  zero-padded  Fourier  tran- 
form.  For  each  of  the  regions  detected  in  the  reference  image,  the  average  magnitude 
of  the  zero-padded  difference  image  is  calculated.  (Note  that  the  number  of  regions 
used  is  not  limited  by  the  number  of  dynamic  encodings.)  The  result  will  be  an  image 
of  the  average  signal  magnitude  change  in  each  of  the  regions.  The  phase  of  the  zero- 
padded  difference  image  is  reintroduced  during  the  creation  of  the  new  constraint 
function  for  the  GS  model.  If  the  original  constraint  function  is  a  baseline  reference 
image,  the  new  reference  image  is 


-f baseline  —  I^baseline  "F  -fdiff,ave  ^ 


(5.7) 


or,  for  the  active  reference  case,  the  new  reference  image  is 

^active  “  -^active  .fdiff,ave  ^ 


(5.8) 


where  /diff,ave  is  the  average  signal  magnitude  change  image  and  Z/diff  is  the  phase  of 
the  zero-padded  difference  image.  This  new  reference  image  is  used  as  the  constraint 
function  in  the  GS  model  to  reconstruct  the  dynamic  image. 


5.2  Results  and  Discussion 

The  proposed  method  was  compared  to  the  original  RIGR  method  with  computer 
simulations  on  experimental  MRI  data  using  both  baseline  and  active  reference  im¬ 
ages.  The  experimental  data  shown  here  are  from  the  contrast-enhanced  dynamic 
imaging  study  of  a  cancerous  rat  described  in  Chapter  4.  As  before,  the  central  N 
encodings  of  the  dynamic  data  sets  were  used  as  the  reduced-encoding  dynamic  data. 
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Figure  5.1:  Explicit  Edge  Information  with  Precontrast  Reference  Image: 
(a)-(b)  The  original  precontrast  reference  and  dynamic  images,  respec¬ 
tively,  reconstructed  using  256  phase  encodings,  (c)-(d)  The  dynamic  im¬ 
age  reconstructed  using  32  dynamic  encodings  with  the  original  RIGR 
method  and  the  proposed  method,  respectively. 

The  results  of  applying  the  proposed  method  with  a  precontrast  reference  image 
are  shown  in  Fig.  5.1.  Images  (a)-(b)  are  the  original  reference  and  dynamic  images, 
respectively,  reconstructed  using  256  phase  encodings.  Images  (c)-(d)  are  the  dynamic 
image  reconstructed  using  32  dynamic  encodings  with  the  original  RIGR  method  and 
the  proposed  method,  respectively.  Note  the  improved  reproduction  of  the  signal 
magnitude  in  the  contrast-enhanced  tumor  with  the  proposed  method. 

The  edge  extraction  step  appears  to  currently  be  the  limiting  step  in  the  process. 
In  the  RIGR  image  of  Fig.  5.1(d),  the  blockiness  of  the  detected  regions  is  apparent, 
in  particular  at  the  outer  edges  of  the  lesion.  This  can  be  reduced  in  some  situations 


52 


Figure  5.2:  Explicit  Edge  Information  with  Postcontrast  Image:  (a)-(b) 
The  original  postcontrast  reference  and  dynamic  images,  respectively,  re¬ 
constructed  using  256  phase  encodings,  (c)-(d)  The  dynamic  image  recon¬ 
structed  with  32  dynamic  encodings  using  the  original  RIGR  method  and 
the  proposed  method,  respectively. 

by  using  an  active  reference  image  in  which  the  contrast  between  adjacent  regions 
is  larger,  making  it  easier  for  the  edge  extraction  step  to  detect  the  edges  well.  In 
addition,  if  there  is  no  edge  information  about  a  new  feature  in  the  reference  image, 
the  edge  extraction  step  will  not  capture  it.  The  result  is  that  the  new  feature 
will  be  averaged  in  with  a  surrounding  area  in  the  contrast  difference  determination 
step.  This  also  supports  the  use  of  an  active  reference  image,  because  it  may  contain 
additional  edges  introduced  by  the  dynamic  changes. 

Both  of  these  effects  can  be  seen  in  Fig.  5.2,  which  shows  the  results  of  using 
edge  information  extracted  from  a  postcontrast  reference  image.  Images  (a)-(b)  are 
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the  original  postcontrast  reference  image  and  the  dynamic  image,  respectively,  re¬ 
constructed  using  256  phase  encodings.  Images  (c)-(d)  show  the  dynamic  image 
reconstructed  with  32  dynamic  encodings  using  the  original  RIGR  method  and  the 
proposed  method,  respectively.  Image  (d)  has  a  sharper  reproduction  of  some  of  the 
internal  details  of  the  tumor  such  as  those  indicated  by  the  arrows.  However,  the 
improvement  between  (c)-(d)  is  not  as  great  as  between  Figs.  5.1(c)-(d)  because  the 
contrast  in  the  postcontrast  reference  image  is  more  similar  to  the  dynamic  image 
than  the  contrast  in  the  precontrast  reference  image. 

5.3  Summary 

A  method  has  been  developed  to  use  explicit  edge  information  extracted  from 
a  high-resolution  reference  image  to  inject  dynamic  information  into  the  GS  basis 
functions.  If  possible,  it  is  preferable  to  use  an  active  reference  image  for  the  edge 
extraction  step,  because  it  may  result  in  better  detection  of  the  image  regions. 
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CHAPTER  6 


MOTION-COMPENSATED  DYNAMIC 

IMAGING 


Object  motion  is  often  a  problem  for  MR  imaging.  For  convenience,  it  is  useful  to 
categorize  object  motion  in  terms  of  its  temporal  and  physical  characteristics.  If  the 
motion  occurs  during  one  Tr  period,  it  is  referred  to  as  intraview  motion;  whereas, 
motion  that  occurs  between  Tr  periods  is  called  interview  motion.  In  addition,  motion 
that  occurs  between  complete  data  sets  is  called  interset  motion.  In  conventional 
MRI,  intraview  and  interview  motion  of  the  object  can  cause  image  artifacts,  because 
it  destroys  the  encoding  relationship  in  the  imaging  equation,  although  the  former  is 
not  usually  a  problem  due  to  the  short  times  involved.  In  addition  to  intraview  and 
interview  motion,  interset  motion  is  also  a  problem  for  constrained  dynamic  imaging, 
because  motion  of  the  object  between  the  reference  and  dynamic  data  sets  can  render 
the  reference  information  useless  as  a  constraint  for  image  reconstruction. 

An  approach  to  overcome  this  problem  is  to  detect  the  object  motion  before  the 
constrained  reconstruction  step  is  performed.  However,  detection  of  object  motion  in 
reduced-encoding  dynamic  imaging  is  nontrivial  due  to  several  factors.  First,  dynamic 
image  contrast  changes  and  object  motion  are  mixed  together.  Second,  the  dynamic 
data  sets  are  low-resolution,  and  it  is  usually  necessary  to  detect  motion  to  a  higher 
accuracy  than  that  dictated  by  the  low-resolution  Fourier  pixel  size.  To  overcome 
these  problems,  we  propose  to  use  a  similarity  norm  which  can  accurately  detect  the 
motion  in  spite  of  the  contrast  changes  and  the  low-resolution  nature  of  the  dynamic 
data.  The  similarity  norm  tries  to  remove  the  effects  of  the  contrast  change  by  using 
only  the  edges  from  the  high-resolution  reference  image  for  the  motion  estimation. 
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6.1  Similarity  Norm-Based  Motion  Estimation 

Given  two  images  Ii  and  I21  we  assume  that 


Ii{xcos6q  -  ysin^o  +  xo,a:sin^o  +  ycos^o  +  Vo)-  (6.1) 


That  is  to  say,  there  is  a  relative  rigid-body  rotation  and  translation  between  Ii  and 
I2,  as  specified  by  xo,  yo,  and  Oq.  The  tilde  signifies  that  Ii  and  I 2  can  have  different 
contrast  behavior.  Therefore,  the  goal  of  the  motion-estimation  step  is  to  find  Xq,  yo, 
and  9o. 

Assuming  that  I2  is  a  low-resolution  image  and  Ii  is  a  high-resolution  image,  we 
first  segment  Ii  into  a  number  of  “homogeneous”  regions.  The  strategy  here  is  to  use 
these  region  boundaries  as  landmark  features.  We  superimpose  this  region  structure 
onto  I2  and  then  calculate  the  regional  intensity  inhomogeneity  af.  We  argue  that 
this  inhomogeneity  is  a  good  indicator  of  the  misalignment  between  the  two  images. 
Specifically,  we  define  the  misalignment  error  Ea  as 


Ea  = 


\ 


Nreg 


i=i 


N 


(6.2) 


where  Areg  is  the  number  of  regions,  mi  is  the  number  of  pixels  in  each  region  and 
N  is  the  total  number  of  pixels.  Clearly,  the  value  of  Ea  is  a  function  of  the  motion 
parameters.  By  minimizing  Ea,  the  values  of  xq,  yo,  and  Oq  can  be  found.  The 
estimated-motion  parameters  that  minimize  Ea  are  then  used  in  the  GS  model  as 


4„{r)  =  T„,/,rt(r)  Y:  (6-3) 

n.eA/jyn 

where  Test  is  the  transformation  that  corresponds  to  the  estimated-motion  parameters. 

The  proposed  method  is  applied  to  a  sequence  of  dynamic  images  in  turn,  so  that 
the  dynamic  data  set  is  compared  to  a  high-resolution  image  that  should  have  a  more 
similar  edge  structure  than  the  original  reference  image.  First,  the  motion  between 
the  high-resolution  reference  image  and  the  first  dynamic  data  set  is  measured.  These 
measurements  are  used  to  correct  and  reconstruct  the  first  GS  dynartiic  image.  This 
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high-resolution  GS  image  is  used  with  the  second  dynamic  data  set  to  determine  the 
motion  that  occurred  between  these  two  acquisition  times.  The  cumulative  motion 
measurements  are  used  with  the  reference  image  to  reconstruct  the  second  GS  dy¬ 
namic  image.  This  procedure  is  repeated  for  the  entire  image  sequence.  Note  that 
the  dynamic  image  is  reconstructed  using  the  original  high-resolution  reference  image 
to  reduce  errors  that  could  arise  due  to  the  accumulated  motion  correction. 

If  it  cannot  be  assumed  that  no  appreciable  motion  occurs  during  the  collection 
of  a  particular  reduced-encoding  d}mamic  data  set,  other  methods  will  have  to  be 
employed  to  measure  this  intraset  motion.  One  possible  way  to  do  that  is  by  using 
navigator  techniques  [91]  that  acquire  a  “navigator  echo”  during  each  view  in  addition 
to  the  image  information.  Each  navigator  echo  is  then  compared  to  an  initial  reference 
navigator  echo  using  correlation  [91]  or  a  least  squares  technique  [92-94]  to  determine 
the  motion  along  the  navigator  direction.  Additional  navigator  echoes  can  be  used 
to  measure  motion  in  the  other  directions,  or  orbital  navigator  echoes  [93]  can  be 
used  to  simultaneously  measure  the  two  directions  of  translation  and  rotation  in  a 
plane.  Although  the  navigator  techniques  are  being  used  in  many  applications,  the 
methods  cannot  be  directly  applied  to  dynamic  imaging,  because  the  basic  assumption 
of  the  navigator  echo  method  will  be  violated;  namely,  that  all  of  the  changes  in  the 
navigator  data  are  due  to  motion  of  the  object.  This  causes  a  problem  in  dynamic 
imaging  in  which  the  navigator  data  from  each  view  can  look  very  different  even 
without  motion,  leading  to  incorrect  motion  estimates. 

One  may  try  to  get  around  this  by  comparing  a  navigator  to  the  immediately 
preceeding  navigator  to  determine  the  motion  parameters  as  opposed  to  an  initial 
reference'  navigator.  Although  this  should  reduce  the  effect  of  the  dynamic  changes 
on  the  navigator  data,  the  incremental  motion  may  be  too  small  to  be  detected  with 
this  method.  Then,  because  the  motion  at  a  point  in  time  would  be  the  accumulation 
of  the  measured  incremental  motions  up  to  that  time,  the  error  in  the  motion  estimate 
may  become  quite  large. 

Perhaps  a  better  way  to  use  the  navigator  method  with  dynamic  imaging  applica¬ 
tions  is  to  design  the  pulse  sequence  such  that  the  navigator  signal  is  not  sensitive  to 
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the  contrast  changes.  For  example,  consider  a  contrast-enhanced  dynamic  study  us¬ 
ing  a  T2  contrast  agent.  The  T2  contrast  agent  modifies  the  appearance  of  the  tissues 
which  take  it  up  by  changing  the  T2  relaxation  constant.  Therefore,  the  image  data 
should  be  sensitive  to  the  change  in  T2,  but  the  effect  of  T2  on  the  navigator  data 
should  be  minimal.  A  possible  way  to  accomplish  this  is  to  acquire  the  FID  signal 
following  the  RF  excitation  pulse  as  the  navigator  data  and  use  the  echo  signal  as  the 
image  data.  In  this  way,  the  image  can  be  T2  weighted,  but  the  navigator  signal  will 
be  proton  density  weighted.  This  method  would  require  careful  design  of  the  pulse 
sequence  with  the  given  application  in  mind. 

In  many  cases,  the  motion  may  occur  in  three  dimensions,  as  opposed  to  the 
planar  motion  discussed  here.  In  this  case,  the  solution  will  depend  upon  whether 
the  imaging  sequence  is  acquiring  2D  slices  or  a  3D  volume.  In  the  case  of  2D  slices, 
the  excitation  and  signal  reception  locations  will  need  to  be  dynamically  adapted 
based  on  the  detected  motion  perpendicular  to  the  imaging  plane  [95].  The  in- plane 
motion  can  then  be  addressed,  as  discussed  previously.  For  3D  imaging,  the  motion 
detection  scheme  would  have  to  be  expanded  to  detect  all  six  degrees  of  motion  (three 
translations  and  three  rotations). 

6.2  Results  and  Discussion 

Computer  simulations  were  performed  to  characterize  the  motion  artifacts  that 
can  arise  with  the  generalized  series  (GS)  model  due  to  interset  motion.  The  sim¬ 
ulations  investigated  the  three  in-plane  motions:  translation  in  the  phase-encoding 
direction,  translation  in  the  frequency-encoding  direction,  and  rotation.  The  ref¬ 
erence  image  shown  in  Fig.  6.1(a)  was  used  for  all  cases  in  the  simulation.  The 
high-resolution  dynamic  image  is  shown  in  Fig.  6.1(b)  for  the  case  of  no  translation 
or  rotation  from  the  reference  position.  Both  of  these  images  were  reconstructed  using 
128  phase-encodings. 

The  first  motion  that  was  investigated  was  translation  in  the  phase-encoding  di¬ 
rection.  Figure  6.2  rows  1-3  show  the  dynamic  image  and  profiles  through  the  upper 
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Figure  6.1:  Motion  Study  High-Resolution  Images:  (a)-(b)  Precontrast 
reference  image  used  for  all  simulations  and  dynamic  image  with  no  trans¬ 
lation  or  rotation  from  the  reference  position,  respectively.  Both  images 
were  reconstructed  using  128  phase-encodings. 


Figure  6.2:  Motion  Study  (Translation  in  Phase-Encoding  Direction): 
(rows  1-3)  Dynamic  images  and  profiles  through  the  upper  and  lower  le¬ 
sions,  respectively,  reconstructed  using  GS  with  32  dynamic  encodings  for 
translations  of  0,  1,  2,  3,  4  and  5  pixels  ((a)-(f),  respectively)  in  the  phase¬ 
encoding  direction  (horizontal)  between  the  reference  and  dynamic  data 
acquisitions. 
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Figure  6.3:  Motion  Study  (Translation  in  Frequency-Encoding  Direction): 
(rows  1-3)  Dynamic  images  and  profiles  through  the  left  and  right  le¬ 
sions,  respectively,  reconstructed  using  GS  with  32  dynamic  encodings 
for  translations  of  0,  1,  2,  3,  4  and  5  pixels  ((a)-(f),  respectively)  in  the 
frequency-encoding  direction  (vertical)  between  the  reference  and  dynamic 
data  acquisitions. 


and  lower  set  of  lesions,  respectively,  reconstructed  using  32  dynamic  encodings  with 
the  GS  model.  Images  (a)-(f)  correspond  to  a  translation  of  0,  1,  2,  3,  4  and  5  pixels, 
respectively,  in  the  phase-encoding  direction  (horizontal)  between  the  reference  and 
the  dynamic  data  sets.  Note  the  ringing  artifact  extending  from  the  right  side  of  the 
phantom.  As  would  be  expected,  the  ringing  is  worse  for  larger  motions  and  for  a 
smaller  number  of  dynamic  encodings.  Note  that  these  artifacts  arise,  even  though 
the  translation  is  smaller  in  all  cases  than  the  low-resolution  Fourier  pixel  size. 

The  simulations  were  repeated  for  translation  occurring  in  the  frequency-encoding 
direction  (vertical)  between  the  acquisition  of  the  reference  and  dynamic  data  sets. 
The  dynamic  image  and  profiles  through  the  left  and  right  set  of  lesions  resulting 
from  a  GS  reconstruction  using  32  dynamic  encodings  are  shown  in  Fig.  6.3  rows  1-3, 
respectively.  Images  (a)-(f)  involve  a  translation  of  0, 1,  2,  3,  4  or  5  pixels,  respectively, 
in  the  frequency-encoding  direction  (vertical)  between  the  reference  and  dynamic 
data  sets.  In  this  case,  the  result  is  also  a  ringing  artifact  emanating  from  one 
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Figure  6.4:  Motion  Study  (Rotation):  (a)-(g)  Dynamic  image  recon¬ 
structed  with  GS  using  32  dynamic  phase-encodings  for  rotations  of  0, 
1,  2,  3,  6,  10  and  90  clockwise  degrees,  respectively,  between  the  reference 
and  dynamic  data  acquisitions. 

side.  However,  the  nature  of  the  ringing  is  slightly  different  from  that  seen  with 
the  motion  in  the  phase-encoding  direction.  This  is  due  to  the  fact  that  the  GS 
model  is  applied  along  the  (low-resolution)  phase-encoding  direction,  and  the  Fourier 
transform  is  applied  along  the  (high-resolution)  frequency-encoding  direction.  As 
before,  the  artifacts  are  worse  for  increasing  motion  and  a  decreasing  number  of 
dynamic  encodings. 

The  last  planar  motion  to  be  investigated  was  rotation.  Figures  6.4(a)-(g)  show 
the  dynamic  image  reconstructed  using  32  phase-encodings  with  the  GS  model  for 
a  rotation  of  0,  1,  2,  3,  6,  10  and  90  clockwise  degrees,  respectively,  between  the 
reference  and  dynamic  data  sets.  In  this  simulation,  a  slight  blurring  in  the  phase¬ 
encoding  direction  of  the  edges  of  internal  structures  results.  The  blurring  increases 
with  increasing  rotation  and  a  decreasing  number  of  dynamic  encodings.  Because  the 
outer  boundary  of  the  phantom  is  symmetric,  no  effect  from  the  rotation  is  seen  in 
the  reconstruction  of  the  outer  boundary. 

To  avoid  these  motion  artifacts,  it  is  necessary  to  detect  and  correct  for  motion 
that  occurs  between  the  acquisition  of  the  reference  and  dynamic  data  sets  prior  to 
constrained  reconstruction.  The  proposed  method  has  performed  well  when  applied 
to  a  sequence  of  dynamic  images,  as  discussed  earlier.  The  results  of  applying  the 
proposed  method  to  a  contrast-enhanced  dynamic  study  of  a  rat  with  cancer  are  shown 
in  Fig.  6.5  in  which  (a)  is  the  reference  image  that  was  reconstructed  using  256  phase- 
encodings.  Images  (b)-(d)  were  reconstructed  using  the  GS  model  with  32  dynamic 
encodings.  In  image  (b),  there  was  no  motion  between  the  reference  and  dynamic 
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(a) 


(b) 


(c) 


Figure  6.5:  Motion  Corrected  Dynamic  Images:  (b)-(d)  were  recon¬ 
structed  using  the  GS  model  with  (a)  as  the  reference  image  and  32 
dynamic  encodings,  (b)  The  dynamic  image  that  results  with  no  mo¬ 
tion  between  the  reference  and  dynamic  data  sets.  The  remaining  images 
show  the  dynamic  image  reconstructed  with  a  5  pixel  shift  in  the  phase¬ 
encoding  direction  (vertical),  a  -3  pixel  shift  in  the  frequency-encoding 
direction  (horizontal)  and  a  3  degree  clockwise  rotation  between  the  refer¬ 
ence  and  dynamic  data  sets,  (c)-(d)  The  dynamic  images  that  result  with 
no  motion  correction  and  with  the  proposed  method,  respectively.  Note 
the  reduced  motion  artifacts  in  (d)  as  compared  to  (c). 
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data  sets,  and  therefore,  it  represents  the  ideal  GS  reconstruction.  The  remaining 
images  represent  a  case  in  which  the  position  of  the  object  during  the  dynamic  data 
acquisition  has  changed  from  the  reference  position  by  a  rotation  of  3  degrees  and 
shifts  of  5  and  -3  pixels  in  the  phase-encoding  and  frequency-encoding  directions, 
respectively.  Figures  6.5(c)-(d)  were  reconstructed  with  no  motion  correction  and 
with  the  proposed  method,  respectively.  Note  the  reduced-motion  artifacts  in  the 
corrected  image  (d). 

6.3  Summary 

The  motion  artifacts  that  are  seen  with  the  GS  model  include  ringing  and  blur¬ 
ring  due  to  in-plane  translation  and  rotation,  respectively,  between  the  reference  and 
dynamic  data  sets.  To  avoid  these  artifacts,  the  motion  between  the  data  sets  needs 
to  be  measured,  which  is  a  difficult  problem  for  many  reduced-encoding  dynamic 
imaging  applications.  The  proposed  solution  uses  a  similarity  norm  to  accurately 
detect  the  motion  in  spite  of  the  contrast  changes  and  the  low-resolution  nature  of 
the  dynamic  data. 
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CHAPTER  7 


FUTURE  WORK  AND  CONCLUSIONS 

7.1  Future  Work 

•  3D  RIGR/TRIGR:  For  many  dynamic  imaging  applications,  it  would  be  de¬ 
sirable  to  do  full  3D  imaging,  as  opposed  to  a  stack  of  2D  slices.  An  example 
is  contrast-enhanced  breast  imaging,  in  which  high-spatial  resolution  in  3  di¬ 
mensions  would  reduce  the  risk  of  a  missed  lesion.  This  would  require  that  the 
RIGR/TRIGR  algorithms  be  extended  to  3D.  In  order  to  do  this,  it  is  neces¬ 
sary  to  solve  a  block-Toeplitz  system  to  determine  the  GS  coefficients.  If  it 
was  desired  to  use  explicit  edge  information  or  motion  compensation,  the  edge 
extraction  technique  and  motion  measurement  method  would  both  need  to  be 
extended  to  3D.  It  may  be  necessary  to  investigate  faster  algorithms  to  perform 
all  of  these  functions  for  the  3D  case  due  to  the  vast  amount  of  data  involved. 

•  Locally  Focused  RIGR/TRIGR:  Given  the  a  priori  knowledge  that  the  dy¬ 
namic  changes  will  occur  in  a  certain  region  of  the  image,  it  may  be  possible 
incorporate  this  information  with  the  GS  model  to  obtain  a  better  reconstruc¬ 
tion  in  that  area  of  the  dynamic  image.  This  would  have  application  in,  for 
example,  following  the  progress  of  a  needle  biopsy  in  interventional  MRI.  This 
work  may  involve  the  development  of  better  methods  of  background  suppression 
to  allow  the  GS  basis  functions  to  only  represent  the  dynamic  changes  in  the 
region  of  interest. 

•  Contrast-Immune  Navigator  Echoes:  As  discussed  in  Chapter  6,  it  may 
be  possible  to  design  pulse  sequences  for  contrast-enhanced  dynamic  imaging 
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sequences  such  that  the  navigator  signal  is  not  affected  by  the  changing  contrast 
due  to,  for  example,  an  injected  contrast  agent. 

•  Adaptive  Data  Acquisition:  If  the  a  priori  information  that  is  available  is 
a  reference  image,  the  use  of  this  information  to  guide  the  data  acquisition  can 
result  in  serious  image  artifacts.  However,  it  may  be  possible  to  use  other  types 
of  a  priori  information  to  guide  the  data  acquisition  in  a  beneficial  way. 

7.2  Conclusions 

The  reduced-encoding  dynamic  MR  imaging  problem  has  been  addressed  with 
the  objective  of  obtaining  simultaneously  high  temporal  and  spatial  resolutions.  To 
avoid  the  loss  of  spatial  resolution  that  occurs  with  reduced-encoding  Fourier  imag¬ 
ing,  many  methods  make  use  of  a  priori  information  at  some  point  in  the  process 
to  reduce  the  number  of  encodings  that  are  necessary  for  the  dynamic  images.  This 
study  reveals  that  if  the  available  a  priori  information  is  a  reference  image,  direct 
use  of  this  information  to  “optimize”  data  acquisition  using  the  existing  wavelet 
transform  or  singular  value  decomposition  schemes  can  undermine  the  capability  to 
detect  new  image  features.  However,  proper  incorporation  of  the  a  priori  information 
in  the  image  reconstruction  step  can  significantly  reduce  the  resolution  loss  asso¬ 
ciated  with  reduced-encoding.  For  Fourier-encoded  data,  we  have  shown  that  the 
generalized-series  (GS)  model  is  an  effective  mathematical  framework  for  carrying 
out  the  constrained  reconstruction  step. 

To  further  improve  the  reconstructed  image,  several  techniques  were  developed 
to  improve  the  GS  basis  functions.  The  two- reference  reduced-encoding  imaging  by 
generalized-series  reconstruction  (TRIGR)  method  uses  a  second  high-resolution  ref¬ 
erence  image  to  suppress  the  background  information.  This  allows  the  GS  basis 
functions  to  more  accurately  represent  the  areas  of  dynamic  change.  A  second  tech¬ 
nique  uses  explicit  edge  information  from  the  reference  image  to  inject  information 
from  the  dynamic  data  into  the  GS  basis  functions.  The  resulting  GS  basis  functions 
more  closely  resemble  those  that  would  be  derived  from  the  dynamic  image  itself. 
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Finally,  motion  that  occurs  between  the  acquisition  of  the  reference  data  set  and 
the  dynamic  data  set  can  render  the  reference  information  useless  as  a  constraint  for 
image  reconstruction.  This  is  a  difficult  problem  to  address,  because  the  dynamic 
changes  may  alter  the  appearance  of  the  image  significantly,  posing  problems  for 
both  navigator-based  techniques  and  registration  algorithms.  A  motion  compensation 
method  is  proposed  which  uses  a  similarity  norm  to  accurately  detect  the  motion,  in 
spite  of  the  contrast  changes  and  low-resolution  nature  of  the  dynamic  data.  The 
technique  was  shown  to  significantly  reduce  motion  artifacts  in  GS  images. 
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APPENDIX  A 


SIGNAL-TO-NOISE  RATIO 


If  the  MRI  measurement  noise  is  assumed  to  be  additive  noise  from  an  ergodic,  sta¬ 
tionary,  white  noise  process  with  zero  mean  and  standard  deviation  ad,  the  resulting 
data  can  be  expressed  as 

d{k)  =  d{k)  +  r]d{k),  (A.l) 

^  * 

where  d{k)  is  the  measured  noisy  data,  d{k)  is  the  noiseless  data  and  r]d{k)  is  the 
measurement  noise.  This  noise  will  be  transformed  to  the  image  domain  during 
image  reconstruction,  resulting  in  the  noisy  image 

J{r)  =  I{r)  +  rii{r).  (A.2) 

Clearly,  the  resulting  noise  image  r]j{r)  will  depend  upon  the  image  reconstruction 
scheme  used,  as  well  as  the  number  of  encodings. 

With  this  in  mind,  the  signal-to-noise  ratio  (SNR)  behavior  of  the  Fourier  series 
and  generalized-series  (GS)  methods  was  investigated  using  ID  profiles.  Both  noise¬ 
less  and  noisy  reconstructions  of  the  dynamic  image  were  created  for  various  noise 
levels  and  various  numbers  of  dynamic  encodings.  The  SNR  of  the  reference  and 
dynamic  data  was  the  same,  reflecting  the  common  experimental  implementation  of 
the  method.  (If  signal  averaging  was  used  to  improve  the  SNR  of  the  reference  data, 
this  would  also  affect  the  SNR  of  the  resulting  GS  dynamic  images.)  The  results  ob¬ 
tained  for  the  reduced-encoding  Fourier  series  and  generalized-series  reconstructions 
were  compared  with  that  of  the  high-resolution  dynamic  image. 

Two  measures  of  SNR  were  calculated.  The  first  measure  included  both  random 
noise  and  systematic  artifacts,  and  the  second  measure  included  only  the  random 
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noise  component.  In  the  first  case,  the  noiseless  dynamic  image  (gold  standard)  was 
subtracted  from,  for  example,  the  noisy  GS  image,  yielding  a  profile  that  included 
both  noise  and  systematic  artifacts.  In  the  second  case,  the  noiseless  GS  image  was 
subtracted  from  the  noisy  GS  image,  leaving  a  pure  noise  profile.  For  each  data 
point,  the  average  standard  deviation  (SD)  of  the  resulting  profile  for  1000  trials  with 
different  noise  realizations  was  calculated  to  quantify  the  performance. 

Figure  A.l  shows  one  realization  of  the  1000  trials  conducted  with  32  dynamic 
encodings.  Rows  1  and  2  show  the  noiseless  and  noisy,  respectively,  reconstructions 
of  the  reference  ((a),(e))  and  dynamic  ((b), (f))  profiles  reconstructed  using  512  en¬ 
codings  and  the  Fourier  series  ((c), (g))  and  GS  ((d), (h))  profiles  reconstructed  using 
32  dynamic  encodings.  Row  3  shows  the  difference  between  the  noisy  reconstruc¬ 
tions  (f)-(h)  and  the  noiseless  dynamic  image  (b).  As  discussed  before,  this  gives  a 
measure  of  the  noise  and  systematic  artifacts  that  are  present  in  the  noisy  reconstruc¬ 
tions.  Row  4  shows  the  difference  between  the  noisy  reconstructions  (f)-(h)  and  the 
corresponding  noiseless  reconstructions  (b)-(d).  Comparison  of  rows  3  and  4  shows 
that  significant  systematic  artifacts  are  present  in  the  Fourier  series  and  GS  profiles 
reconstructed  with  32  dynamic  encodings  in  this  simulation.  Note  the  reduced  noise 
level  in  the  truncated  Fourier  reconstruction  (m)  due  to  the  reduced  number  of  data 
points  that  are  used  in  the  reconstruction. 

Figure  A. 2  shows  the  mean  SD  of  the  noise  alone  for  the  Fourier  series  and  GS 
reconstructions  as  the  number  of  dynamic  encodings  ranges  from  8  to  256  of  512 
encodings.  As  expected,  the  Fourier  series  reconstructions  have  a  lower  mean  SD  at 
every  number  of  dynamic  encodings  than  the  high-resolution  dynamic  reconstruction 
due  to  the  reduced  number  of  data  points  used.  The  mean  SD  for  the  GS  model  is 
consistently  higher  than  that  of  the  high-resolution  dynamic  reconstruction.  For  this 
simulation,  the  mean  SD  of  the  GS  reconstruction  improves  only  slightly  as  a  greater 
number  of  dynamic  encodings  is  used  past  32  encodings. 

Figure  A. 3  depicts  the  behavior  of  the  mean  SD  of  the  noise  and  systematic 
artifacts  for  the  Fourier  series  and  GS  reconstructions  with  a  varying  number  of  dy¬ 
namic  encodings  ranging  from  8  to  256  of  512  encodings.  The  mean  SD  for  the 
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Figure  A.l;  SNR  Study:  Rows  1  and  2  are  the  noiseless  and  noisy  recon¬ 
structions,  respectively,  of  the  reference  ((a),(e))  and  dynamic  ((b), (f)) 
profiles  reconstructed  with  512  encodings  and  the  truncated  Fourier 
((c),(g))  and  GS  ((d),(h))  profiles  reconstructed  using  32  dynamic  en¬ 
codings.  Row  3  shows  the  difference  between  the  noisy  reconstructions 
(f)-(h)  and  the  noiseless  dynamic  profile  (b).  This  gives  a  measure  of  the 
noise  and  systematic  artifacts  that  are  in  the  noisy  reconstructions.  Row 
4  shows  the  difference  between  the  noisy  reconstructions  (f)-(h)  and  the 
corresponding  noiseless  reconstructions  (b)-(d).  This  gives  a  measure  of 
the  noise  performance.  Note  that  the  four  rows  are  not  on  the  same  scale. 
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Mean  of  Standard  Deviation  of  1000  Trials 

Difference  with  Corresponding  Noiseless  Image 


0.0  100.0  200.0  300.0 
Number  of  Dynamic  Encodings 

Figure  A. 2:  SNR  Behavior  of  Noise:  Curves  show  the  dependence  of  the 
mean  SD  of  the  noise,  with  respect  to  the  number  of  dynamic  encodings 
used  in  the  reconstruction  for  the  Fourier  series  and  generalized-series. 
The  mean  SD  from  the  high-resolution  dynamic  image  is  included  as  a 
horizontal  line  for  comparison  purposes. 
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Figure  A. 3:  SNR  Behavior  of  Noise  and  Systematic  Artifacts;  Curves  show 
the  dependence  of  the  mean  SD  of  the  noise  and  systematic  artifacts,  with 
respect  to  the  number  of  dynamic  encodings  used  in  the  reconstruction 
for  the  Fourier  series  and  generalized-series.  The  mean  SD  from  the  high- 
resolution  dynamic  image  is  included  as  a  horizontal  line  for  comparison 
purposes. 
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high-resolution  dynamic  image  is  included  as  a  horizontal  line  for  comparison  pur¬ 
poses.  As  would  be  expected,  the  combination  of  noise  and  systematic  artifacts 
increases  as  the  number  of  dynamic  encodings  decreases  for  both  the  Fourier  series 
and  the  generalized  series.  The  GS  reconstructions  have  a  smaller  mean  SD  than  the 
Fourier  series  reconstructions  at  all  numbers  of  dynamic  encodings  due  to  the  reduced 
systematic  artifacts,  as  can  be  determined  using  both  Figs.  A. 2  and  A. 3. 

In  summary,  the  noise  behavior  of  the  reduced-encoding  Fourier  series  is  better 
than  that  of  the  generalized-series.  However,  if  both  random  noise  and  systematic 
artifacts  are  considered,  the  truncated  Fourier  series  no  longer  performs  better  than 
the  generalized-series.  Because  the  quality  of  the  resulting  dynamic  image  is  affected 
by  both  random  noise  and  systematic  artifacts,  the  generalized-series  would  be  the 
better  choice  for  reduced-encoding  dynamic  image  reconstruction. 
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